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ABSTRACT 



A binary collision simulation of the radiation damage problem is 
proposed. Energy thresholds and parameter adjustment are discussed. 
Preliminary results show strong focusing in the (110) direction and 
weaker focusing in the (100) direction. 

Although not complete, the model is expected, with suggested 
modifications, to adequately represent the radiation damage problem. 
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and encouragement given them by Associate Professor Don E. Harrison, Jr. 
of the U. S. Naval Postgraduate School in this investigation. The authors 
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1 . Introduction 



The interaction of energetic particles with matter may be divided 
into two energetic regions of interest: 

1) High energy events in which interactions are primarily 
inelastic electronic processes of excitation and ionization, 
and 

2) Low energy events in which the interactions are essentially 
elastic interactions between whole atoms. 

There is an intermediate energy range in which both processes are im- 
portant, but until the two limiting cases are better understood, attempts 

1 * 

to understand the intermediate range seem futile. The threshold energy 
for which electronic excitation and ionization processes become impor- 
tant depends upon the mass of the energetic particle. An approximate 
formula for this threshold energy is E = A kev. , where A is the atomic 
mass of the energetic particle. The threshold for electronic processes 

seems to be relatively insensitive to the atomic mass of the target 
2 

material. For metals this threshold is of the order of 10-100 kev. The 
high energy processes of electronic excitation and ionization have re- 
ceived considerable attention and are reasonably well understood. 

We will confine our attention to the energy region below this threshold 
for inelastic electronic processes. 

Many physical properties of a material are changed when the 
7 

material is irradiated. The Young's modulus and the electrical 

* All footnotes refer to the bibliography. 
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resistivity are among those properties of materials which show most 
significant changes., For metals and semiconductors, the residual 
electrical resistivity is most sensitive to radiation. 

Semiconductors, for example, may show changes in electrical 

7 

resistivity of up to four orders of magnitude, but there is no known 
simple relation between the magnitude of the change in electrical 
resistivity and the amount of radiative energy absorbed by a semi- 
conductor. Radiation changes some semiconductors from n-type to the 
p-type . 

The radiation damage problem involves chemical bonding between 
the atoms of the target material as well as elastic collisions between 
atoms. For low energy collisions, chemical bonding forces vary widely 
for different atomic species, but are of the same order of magnitude as 
the elastic repulsive forces. Therefore, to reduce the complexity of 
the problem, we have limited our investigation to simple crystals of 
face centered cubic metals, with special emphasis on copper. 

Davies, et al. have developed a chemical technique for removing 

g 

thin layers (2> 37 A) from aluminum foils. The thickness of the layers 
is easily controlled. The technique has been widely used to study the 
range of radioactive ions in aluminum.^ Ranges have been measured 
for Cs 137, Rb^S, Na^, and K^. For example, the mean range of 
10.5 kev. Na^4 ions in aluminum was found to be about 250 Angstroms. 

When a metal is irradiated with heavy energetic ions ( > 50 kev.) a 
significant number of atoms are ejected from the target through the surface 



2 



which is bombarded by the incident ion beam. This phenomenon is 
known as sputtering. Under suitable conditions as many as 40 atoms or 
ions may be sputtered for each incident particle. ^ Examination of the 
angular distribution of sputtered particles reveals definite preferential 

ejection in directions which correspond to (110), (100), and (ill) lattice 

, 12-16 
directions . 
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2. Lattice Effects 



Whether a lattice atom leaves its original site or not after a 
collision depends upon the energy absorbed by the atom and its final 
direction of motion with respect to the lattice. If an energetic lattice 
atom does not leave its site, it will transfer energy to its neighbors 
through elastic collisions until thermal equilibrium is reached. Be- 
cause of crystalline structure, the dissipation of energy is expected to 
be anisotropic. The only exception to the eventual thermalization of 
the excess energy of a bound lattice atom is the case in which sufficient 
energy reaches a surface of the material to knock one or more atoms away 
from the surface . 

If a lattice atom is given sufficient energy to leave its site per- 
manently, and no other atom replaces the energetic atom, a vacant 
lattice site is created. The surrounding lattice atoms relax slightly 

toward the vacant site, but the lattice structure remains essentially the 

17 

same as that of a perfect lattice. The energetic atom travels through 

the lattice giving up energy through elastic collisions with the atoms 

near its path until it no longer has sufficient energy to force its way 

between the atoms of the crystal. If there are no vacant lattice sites 

near the position at which an energetic atom is brought to rest, it will 

be forced to occupy some position other than a regular lattice site. 

Three possible equilibrium positions have been suggested for these 

interstitial atoms in a face centered cubic crystal: the crowdion , the 

17 

split interstitial pair, and the center of a unit cell. The crowdion 
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consists of an extra atom in a (110) line of atoms. In the split inter- 
stitial pair configuration, the interstitial shares a regular lattice site 
with another lattice atom so that each is symmetrically displaced from 
the equilibrium position in a (100) direction. An atom in the center of 
a unit cell may also be described as a crowdion in a (100) direction. 

The lattice atoms surrounding an interstitial are displaced from their 
regular lattice sites to accommodate the extra atom. More will be said 
about the interstitial configuration later. 

An energetic atom (knock-on) may transfer essentially all of its 

18 

energy to a single lattice atom in a head-on collision. In this case, 
the lattice atom becomes the knock-on, leaving a site for the former 
knock-on to occupy. The new knock-on may, in turn, transfer essentially 
all of its energy to another atom. Thus, a chain of replacements may 
transport an interstitial to a position which is several times the dis- 
tance between nearest neighbors from the corresponding vacancy. Re- 
placement chains may be cyclic so that the last replacement in the chain 

17 

fills the initial vacancy, leaving no vacancies or interstitials. The 
configuration consisting of a vacancy and a nearby interstitial, known 
as a Frenkel pair, may be stable or unstable. The stability of a Frenkel 
pair depends upon both the distance of separation and the orientation 

(with respect to the lattice) of the interstitial and the vacancy. 

19 

Silsbee has shown that an isolated line of hard spheres may 
focus and transfer energy along the line. The condition for focusing 
is that s/d<. cos ^ , where s is the distance between adjacent spheres, 
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d is the diameter of the spheres , and @ is the angle between the 
velocity of the first atom in the line and the line. The hard core 
approximation (diameter of spheres equal to distance of closest 
approach for a head-on collision) is useful in this work. For reason- 
able repulsive potentials such as the Coulomb, Bohr, Born-Mayer, etc. , 
this modified model predicts a focusing effect below a threshold energy, 
and a de-focusing effect above this threshold energy.^ When the 
Gibson potential #2 simulates copper, this model predicts focusing in 
the (110) directions below about 60 ev. , with sharp focusing below 
about 15 ev. The presence of adjacent lines of atoms will enhance 
the focusing effect found in an isolated line of atoms . Although the 
hard sphere models used to demonstrate the phenomenon of energy 
focusing contain rather drastic assumptions, the same sort of effect is 

expected from models based on more realistic potentials. 

20 

A channeling effect has also been reported. This is the tendency 
for an energetic atom traveling in a (110) or (100) direction to be focused 
along a (110) or (100) line that is centered between surrounding (110) or 
(100) lines of lattice atoms. The energetic atom may travel a large dis- 
tance along a channel compared to the distance it would travel in some 
other direction. The rate of energy loss along the channel is highly 
dependent on the potential function which governs the interaction , the 
energy of the atom, and its deviation from the line of focusing. From 
geometrical considerations, one would expect (100) channels to dissipate 
more energy per unit length than would be dissipated by (110) channels. 
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Channels, energy chains, and the displacement chains, then, may 
produce vacancies which are separated by large distances from their 
corresponding interstitials and sputtered atoms. Thus, chains and 
channels are competing mechanisms which produce the same effects . 
Energy and displacement chains should not be significantly affected by 
the presence of vacancies. Vacancies in the lines of atoms surrounding 
a channel should tend to divert channeled energetic atoms into the 
corresponding chains. Interstitial atoms should tend to disrupt channels 
and chains . 

It appears that the chain and channel mechanisms are the dominant 
factors in the radiation damage event. The conditions for stability and 
the lengths of the chain/channels will then largely determine the final 
lattice configuration. 
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3. Theoretical Work by Others 

17 

Gibson, et al„ have studied radiation damage events by inte- 
gration of the classical equations of motion of the atoms in a small 
crystallite (<C1000 atoms) . A central repulsive potential of the Born- 
Mayer type, A exp (-r/b), was assumed for the interactions between 
atoms. In order to contain the events it was necessary to restrict the 
energy to <400 ev. Calculations were made on a high speed digital 
computer. 

17 

The computer program developed by Gibson, et al. keeps a 
record of position co-ordinates and velocity components for each atom 
in the crystallite. The orbits of the atoms are arbitrarily divided into 
segments with respect to time, i.e, time steps. In order to calculate 
the trajectories of the atoms for a particular time step, the atoms are 
considered sequentially. The resultant force on an atom under con- 
sideration is calculated from the position of the atom and the positions 
of the surrounding atoms . The atom is then moved for one time step 
under the influence of this resultant force. When the calculations are 
completed for all atoms for a given time step, calculations for the next 
time step are begun. 

If the repulsive force between adjacent atoms were the only force 
simulated, the crystallite would expand without limit. Therefore, it is 
necessary that an attractive force act on the surface atoms of the crystal- 
lite. If all collisions were elastic collisions, and if the original energy 
(for a 100 ev. run) were completely contained in the lattice, the mean 
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energy of the atoms would be about 0.1 ev. after thermalization . 

This corresponds to about 900° C. as the final temperature of a sys- 
tem which started at 0° C. Therefore, a viscous damping force was 
applied to surface atoms to dissipate energy. 

This model is useful in studying low and moderate energy (<400 ev.) 
effects in radiation damage events . The effects studied were stability 
of lattice defects , anisotropicity of energy transfer (energy chains) , 
and displacement chains. Interatomic potential and surface effects 
were also studied. 

In the Gibson model, energy chains play a dominant role in the 
dissipation of energy from the initial knock-on. Frenkel pairs are 
shown to be stable for distances of separation greater than 2 1/2-4 
times the nearest neighbor distance, depending upon the orientation 
of the pair with respect to the lattice. Interstitials were studied ex- 
tensively and found to be stable only in the split interstitial configura- 
tion . 

21 

Oen, et al. have investigated the stopping of energetic atoms 
(1-100 kev.) by solids. The model is based on binary elastic collisions 
between a primary knock-on and successive target atoms which are 
initially at rest. Lattice effects are neglected; the target atoms are 
initially randomly distributed throughout the target material. 

The ultimate goal of the work by Oen, et al. is to discover the 
appropriate interatomic potential. The immediate aim is to calculate 
histories of many primaries which may be considered to simulate a 
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beam of incident atoms or ions. 



Trajectories and histories are calculated only for the primary 
knock-ons , i.e. any effect produced by sedondary knock-ons is 
neglected. The history of each primary consists of: 

1) R, the distance from the initial position to the final 
position of the primary, 

2) X, the component of R parallel to the initial direction of 
motion of the primary, 

3) P, the component of R in a plane perpendicular to the 
initial direction of motion of the primary, and 

4) L, the length of the path traversed by the primary. 

Each run consists of the calculation of histories for enough primaries, 

usually 1000, to insure that statistical fluctuations in the mean values 

of R, X, P, and L were acceptable (< 3%). 

The potential used was an exponentially screened Coulomb, or 

Bohr potential: V<r) = e *P 

* & b 

where and a b= K 4 h /( 2 ,$ -1“ 7 ^^ 

where a^ = first Bohr radius (0.529 A) . The parameter k may be adjusted 
to achieve agreement with experimental results. The hard sphere approxi- 
mation to the Bohr potential was also used in some cases for comparison. 
Results for the cases of 10.5 and 60 kev. Na^ ions incident upon 

aluminum targets were compared with the experimental results of Davies & 

0 

Sims. No serious discrepancy is found in the mean penetration , X . 
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However, the experimental distribution of penetration is more "skewed” 
toward lower penetration than the distribution of penetration calculated 
by Oen, et al. 
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4. The Problem 



The binary collision model has been the basis of all analytical 
studies of radiation damage and sputtering. Many machine calculations 
have been made in recent years using an N - body collision model. 
However, no serious attempt has been made to determine whether or 
not a two-body approximation would yield satisfactory results . 

Our problem, then, was to build a model based on the binary col- 
lision whose logic would be sophisticated enough to adequately repre- 
sent the radiation damage problem. The nature of the results obtained 
would then definitely establish the degree of confidence to be placed 
in any binary collision approximation models. 
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5. The Binary Collision 



With minor modifications, any repulsive potential may be used in 
the two body interaction. An eroded form of the potential is used, i.e. 
the potential has a definite radius of effect. We have taken this radius 
of effect (or sphere of influence) as an input parameter, SPHI . There is, 
then, no interaction between atoms whose separation is greater than 
SPHI. 

All applications of the program to date have used an eroded Born- 
Mayer potential, V = A exp (B (l — r /SPHI)) for r< SPHI. The constants 
A & B were taken from Gibson et al. for three different Born-Mayer 
potentials, termed simply Potentials, 1,2, and 3. The constants for 
these three potentials are: 

POTENTIAL A (ev.) B 

1 .0392 16.97 

2 .0510 13.00 

3 .1004 10.34 

The interaction starts with the two atoms at SPHI distance of sepa- 
ration with all the energy as kinetic energy. All co-ordinates and velo- 
city components are transformed to the center of mass system, where 
the equations of motion are: 

1) dt=(2/m' (E-V-L 2 /2m'r 2 )) -1 / 2 dr 

2) d0 = (L/m'r 2 ) dt 
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where t = time 

m' = reduced mass 
E = kinetic energy 
V = potential energy 
L = total angular momentum 
r = distance of separation 
0 = angle of deviation 
s = impact parameter 

The equations are integrated from r = SPH1 to r = CPA, the distance 
of closest approach. Since the interaction is symmetrical in the center 
of mass system, the definite integrals of the above equations are one- 
half the time of interaction and one-half the total angular deviation, 
respectively . 

22 

CPA is found by an interation of the formula: 

(CPA) 2 (1 - V(CPA/E r ) - S 2 = 0 

The first approximation to CPA is taken as SPHI. A four point Gauss- 

23 

Legendre quadrature is employed to evaluate the integrals. 
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6 . Development of the Model 



When the possibility of a feasible two body approximation to the 
radiation damage problem is considered, two major problems are im- 
mediately encountered . 

The atoms would, of necessity, be considered in some sort of 
ordered sequence, but all calculations for a single atom could not be 
carried through before the other atoms are considered since this would 
obviously leave out of the calculations the effect of collisions suffered 
by other atoms. This first major problem led to the time step approach, 
and eventually to the use of an absolute time of an atom's last collision. 

The second major problem is to develop a mechanism which will 
handle simultaneous collisions within a two-body approximation. 

After the program was developed, we discovered that the lattice 
would not hold together which necessitated a further modification. 

This difficulty is only partially resolved. 

The first question considered in the development of the model was 
to mathematically create a face centered cubic lattice in table form, 
together with a method which tags each atom so that its subsequent 
movements and initial position could be found. Figure 2, APPENDIX II, 
shows the ordered arrangement found in a (100) plane of a fee lattice. 

The distance between two atoms along the (100) direction is taken as 
two units. It will be noted that if the total X dimension, IX, is odd, 
then there are an equal number of atoms in the line y = 0 and the line 
y = 1 . By specifying EX odd , the number of atoms in a line y = constant 
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is known to be (EX + l)/2. This value, NLINE, can then be multiplied 
by the total Y dimension plus one (IY+1) , and the number of atoms in 
the whole plane, NPLANE, is obtained. It will be noted from Figure 3 
that the next plane, one Z unit removed, is akin to a photographic 
negative of the plane under consideration. However, the ideas of 
NLINE and NPLANE are still valid . The total number of atoms in the 
lattice can then be found by multiplying NPLANE by the Z dimension 
plus one (IZ+1). If the atoms are numbered sequentially as shown in 
APPENDIX II, the initial co-ordinates of their sites can also be found 
by the method shown in APPENDIX II and APPENDICES IV & V (BOX 31) . 

The process outlined above will create and number a 1500 atom 
lattice in approximately one second. An advantage of this method is 
that a rectangular lattice of any size or shape desired can be created 
by specifying IX, IY, and IZ. The only restriction is that IX must be 
odd . 

The atoms are now considered sequentially, for a period of one 
Time Step Length (TSL) , and all collisions that will take place before 
the end of the time step are allowed to happen . All atoms are then 
moved to points in space and time corresponding to the end of the time 
step, and the process is repeated for the next time step. 

The length of a time step was originally calculated so that the 
most energetic atom in the lattice could move one lattice unit (1.807 
Angstroms for copper) in one time step. This has been modified so 
that it is now able to move any desired fraction, TFAC , of one lattice 
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unit. TFAC is, naturally, one of the input parameters of the problem. 

Since the energy of the most energetic particle, ENERGY, in the 
lattice decreases during the course of the run, the TSL should be re- 
calculated periodically for greater efficiency. Provisions were made 
to recalculate the Time Step Length every JPB (an input parameter) time 
steps. The TSL is also recalculated every time another particle is shot 
into the lattice, since it is assumed that this particle will be the most 
energetic one present, and it was not thought feasible to wait for the 
next regular recalculation. A table of the first 30 Time Step Lengths, 
and the time step of recalculation, is kept for reference. 

As now constituted, there are one external and two internal methods 
used to determine the number of time steps which complete a run. The 
external method is the only one available to the operator once a com- 
putational run has commenced. 

First, a variable limit is set on the maximum number of time steps 
(MNTS , an input parameter) that may be taken . The calculation must 
cease when this number of time steps has been completed, and the 
results are then obtained . 

Second, the program, at the beginning of each time step, decides 
whether or not it should continue for another time step. The decision 
to continue is made if the energy of the most energetic particle (called 
ENERGY) in the lattice (excluding all those that have left the lattice) 
is greater than a certain threshold energy, THERMAL (an input parameter 
explained below) . If ENERGY is lesli than THERMAL, the run is 
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considered completed and the informational output is given. 

Third , the operators may decide to shut down operation before 
either of the two internal limits have been reached. This may be done 
at any time, simply by throwing a switch, and the computation will 
cease at the beginning of the next time step. 

The external mode of problem completion should not be the normal 
operating condition, nor should the number of time steps be allowed to 
approach MNTS under most conditions, since either method still leaves 
atoms in the lattice with more energy than the quantity THERMAL. 
Usually the calculation should cease of its own accord, with the 
second method given above as the basis for the internal decision. 

Other factors, such as running time, mechanical failures, etc. , may 
necessitate the use of the external method, but these are to be avoided 
if possible . 

A method of program regeneration was also devised so that the run 
could be discontinued at the start of any time step and then restarted at 
that point at some later time. This has definite value where continuous 
periods of available machine time are not long, or when mechanical or 
electrical failures destroy the calculations . The information needed for 
regeneration can be obtained at the beginning of any time step and the 
calculation continued (or stopped if desired) . 

In this manner, if the calculation has been running for three hours 
before a failure of some kind occurs , the only time that need be re- 
peated on the machine is the time elapsed since the last regeneration 
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information was put on magnetic tape. It should be noted that the re- 



generation need not be undertaken at once , but may be delayed until 
convenient. This regeneration process has great practical importance 
to the user of the program, but is of no theoretical value. 

Many atoms are given only a cursory examination during a time 
step, while the more energetic atoms are subjected to a highly detailed 
treatment. The use of two energy thresholds in the model simplifies 
this choice for a particular atom. 

ETHRESH, the upper energy threshold, is the energy required for 

an atom to move from its site to a site near it, i.e. the displacement 

energy. Calculations by others ^ have shown this displacement 

threshold to exist, and also that it is in the neighborhood of 20 

1 7 

electron volts. They have also shown that this displacement thresh- 
old depends upon direction, and is lower in the (100) and (110) direc- 
tions than in other directions in the lattice. This directional depend- 
ence of the displacement threshold has not been built into the present 
model, but it should be considered in any future development of the 
model . 

THERMAL, the lower energy threshold, is an arbitrary dividing 
line imposed upon the model, below which atoms are quite arbitrarily 
placed on sites or formed into interstitial pairs with other atoms and 
not allowed to move. Atoms above THERMAL are allowed to participate 
normally in collisions . 
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There are, then, three distinct energy classes among the atoms in 



the lattice: 

1) Those with energy above ETHRESH, 

2) Those with energy between ETHRESH and THERMAL, and 

3) Those with energy below THERMAL. 

The three classes will be discussed separately. 

The atoms above ETHRESH do not occupy a site in the lattice, 
collide with other atoms as necessary, and obtain information about 
their subsequent motion only from the interaction subroutine . 

The atoms between ETHRESH and THERMAL are forced to occupy 
specific lattice sites, but are allowed to vibrate around these sites. 
They cannot move to a new site because their energy is less than the 
displacement energy. These atoms are allowed to move undisturbed 
until they have a collision, and then after the collision has been com- 
pleted, if they lose a specific percentage, TURN (an input parameter), 
of their original energy their velocity vectors are changed to point back 
toward the site occupied by the atom. 

We found prior to the use of this turn-around method that the lattice 
had no inherent cohesiveness. This method holds the lattice together. 

To speed up computation time, atoms with energy below THERMAL 
are forced to occupy a fixed site, or to share an interstitial site with a 
partner, and are not allowed to move. This is a reasonable and valuable 
approximation for atoms of low energy, since they will not have an 
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appreciable effect on the results of the computation. The energy and 
associated velocity components of these thermalized atoms are kept 
intact for use in possible future collisions. 

After the decision is made to continue for another time step, and the 
Time Step Length has been recalculated, the atoms are considered sequen- 
tially, The information associated with the atom under consideration 
(primary atom) which is stored in the LB array is first inspected. If the 
atom has: 

* • • i 

l 

1) An energy equal to or less than THERMAL, or has 

2) Been through the current time step , or has 

3) Left the lattice , 

we move on to consideration of the next atom. 

We will account for all collisions if we consider only those atoms 
with energy greater than THERMAL, since atoms with energy less than 
THERMAL are not allowed to move unless hit. Atoms that have already 
been through the time step have completed all interactions which occur 
during the time step and should not become involved again until the 
next time step. Atoms that have left the lattice will not be considered 
for the remainder of the problem . 

A primary atom with energy above THERMAL is finally selected, or 
the computation shuts down. A list of all other atoms it is possible for 
this atom to hit in this time step is required. To shorten the computa- 
tion , a mathematical box is constructed around the atom in question 
(hereafter referred to as atom M) and a list (the NUM list) is made of 
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all atoms whose positions (centers) are inside the box. If the box is 
made large enough, it will be impossible for atom M to hit any atoms 
outside the box, so we need consider only those on the NUM list. At 
present, the box is given sides of 3.02 units, orjf 1.51 units. Initially, 
the sides were + 1.01, but this has been found to be too small, and 
obvious collisions were missed. Atoms that have been through the time 
step are not included on this list, since their co-ordinates are relative 
to the end of the time step rather than the start of the time step. Atom M 
is placed first on the list for programming facility and easy reference. 

The impact parameter between every other atom on the NUM list and 
atom M is then calculated, these distances comprise the DSTANCE list. 
The time at which each atom reaches this distance of closest approach, 
and the time at which the atom and atom M are SPHI (sphere of influence 
of an atom) distance apart are also calculated. These times are placed 
on the TMIN and T lists, respectively. 

During the course of the calculation of the DSTANCE, TMIN, and T 
lists, several tests are performed on these quantities. If the DSTANCE 
for an atom is larger than SPHI, no collision will occur. If the relative 
velocity is less than 10“'®, then the atoms are moving very slowly 
relative to one another, and the collision is ignored. Several tests are 
then performed on each T. If the time T is less than the start of the pre- 
ceeding time step, it would have been considered in this previous time 
step, and is therefore treated as a spurious collision. If the time T is 
greater than the end of the current time step, it will be re-examined in 
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the next time step. The times of the last collision of atom M and the 
other atom are then checked, and if the time T, on an absolute scale, 
is less than either of these, the collision is spurious since it would be 
starting before one that had already taken place. The last test performed 
is whether the time T is within .00001 (TSL) of the end of the time step. 
If so, we feel that it would be better to consider this collision in the 
next time step, since to consider it now would mean that the end of the 
collision was far into the next time step, which is unsatisfactory from 
a practical viewpoint. 

If, after this series of eliminations, there are no atoms left, then 
atom M is moved to the end of the time step, and the next atom in se- 
quence becomes atom M. 

If atoms remain on the list we select the ones that atom M will hit. 
Two alternative methods of selection were considered by the authors: 

1) Select the one, or ones, that are closest (i.e. a space- 
wise selection) , or 

2) Select those that will collide with atom M first (i.e. a 
time-wise selection) . 

Only the second method has been used. 

The space-wise selection would surely include all the hardest 
collisions, and would probably be easier from a programming viewpoint. 
However, a time-wise selection insures that atom M will hit first what 
it really should hit first. Atom M may thereby miss a hard collision 
with another atom because of another previous soft one , but we 
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currently feel that this is a truer picture of what actually happens in the 
lattice during the radiation damage process. We feel that both methods 
should be tried in future development of the model. A study as to the 
relative feasibility of both methods would enhance the usefulness of the 
model . 

Atoms are selected that have the smallest times , T. Since we assume 
that two atoms have no interaction if they are separated by a distance 
greater than SPHI, the time T is the actual time of collision. We feel 
that some leeway should be available here , therefore all atoms whose 
times T are within a certain amount (CUTOFF) of a time step from the 
minimum times are also selected. These atoms are placed on the KHIT 
list, with atom M first on the list. Since CUTOFF is an input parameter, 
adjustment to a satisfactory value may be accomplished easily. Use of 
the CUTOFF parameter also eliminates some of the objections to the 
time-wise selection process. 

The initial velocity components and energies of all atoms on the 
KHIT list are then stored in the SAVE array because they will be needed 
for the velocity scaling processes. 

It must be emphasized that atom M hits all those on the KHIT list 
almost simultaneously „ This necessitates some sort of energy and 
velocity scaling process, or a general solution to the N body problem. 

We feel that if the two body approximation is to have any validity 
atom M must enter into each of the two body collisions with its initial 
energy, but the program must scale the results in some fashion. Assume, 
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for example, that M initially has an energy of 100 ev. and hits two 
atoms simultaneously with identical impact parameters. It is more 
feasible to assume that M hits each one with an energy of 100 ev. 
than to assume that M hits one with 100 ev. and hits the other with 
whatever energy is left, say 37 ev. 

Naturally, with M hitting each atom separately with 100 ev. the 
struck atoms will receive too much energy, we must therefore scale the 
results . 

In the actual collision process of the model, then, atom M and the 
next atom on the KHIT list are moved to their calculated positions at the 
time T for the atom being hit. The collision is then calculated and the 
struck atom is given its new position co-ordinates, velocity components, 
and energy. Atom M and the next atom on the KHIT list are then moved 
to the positions corresponding to their actual time of collision, T, and so 
forth. This is repeated until atom M has "struck" all the atoms on the 
KHIT list. 

To conserve both energy and momentum in any scaling process 
would constitute a general solution to the N body problem. We, there- 
fore, elected to conserve energy rather than momentum, feeling that 
this would give us a more realistic insight into the radiation damage 
problem. 

There are four general scaling cases that must now be considered: 

1) Atom M "lost" more energy than it started with, 

2) Atom M "lost" exactly all the energy it started with, 
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3) Atom M started with an energy above ETHRESH and did not 

"lose" all of its energy, and 

4) Atom M started with an energy below ETHRESH, but above 
THERMAL, and did not lose ail of its energy. 

These cases will be discussed in order. 

CASE I 

Atom M is allowed to lose all its initial energy, and its velocity 
components and energy are set equal to zero. The energies and 
velocities of the struck atoms are then scaled to conserve energy. The 
energy scaling factor used in this case is simply the total energy (of 
those on the KHIT list) before the collisions divided by the total energy 
after the collisions. The velocity scaling factor is the square root of 
the energy scaling factor. The energies and velocities of all atoms on 
the KHIT list , except atom M , are scaled by these two scaling factors . 

Both energy and the absolute velocity directions are conserved by 
this process. 

CASE II 

Atom M is allowed to lose all of its initial energy, and its velocity 
components and energy are set equal to zero Here, however, energy 
scaling is not necessary since energy has already been conserved. 

CASE III 

Here it is not necessary to scale the energies of the struck atoms, 
since atom M did not lose all of its initial energy, so the struck atoms 
retain the energies and velocities they received in the two body collision 
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process. The energy of M is found by subtracting the energy gained by 
the struck atoms from M's initial energy. M’s velocity components 
must now be found , 

The sums of the individual velocity components calculated for M, 
but not given to M, are taken. Call these sums S^, S2, and S3 for the X, 
Y, and Z directions. Let R be (Sj + The magnitude of M's 

final velocity, V, is known since the energy is known. The X velocity 

I 

is then taken as (V/R) (Sj), and similarly the Y and Z velocities are 
(V/R)^) and (V/RHS3) respectively. 

CASE IV 

Two sub-cases must be considered. First, if atom M does not lose 
TURN % of its initial energy, the treatment used in CASE III will be 
followed . 

Second, if atom M does lose a significant portion of its energy, its 
velocity components are changed so that the atom is directed back toward 
the site it occupies . The energy of M is again its initial energy minus 
the energy gained by the struck atoms. The scale factor used is 
(Velocity/distance from site). This is multiplied by AX, AY, and 
A Z to give the velocity components . ^X is simply the X co-ordinate 
of the site minus the X co-ordinate of M. 

In all cases, atom M's spatial co-ordinates are found by averaging 
the individual co-ordinates calculated for M in the X, Y, and Z directions 
(i.e. dividing the sum by NHIT) . 
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If atom M is a primary atom (i.e. one of those considered sequen- 
tially for this time step) , then the KHIT list is duplicated by the LAST 
list. Atom M will then proceed to go through its designated collisions, 
look for more collisions, and go through these until it completes all 
collisions in this time step. 

When atom M has completed all its possible collisions for this time 
step, each atom that M hit (found in order on the LAST list) undergoes 
the same process and attempts to hit other atoms. The atoms hit by M, 
or others on the LAST list, are not inspected for possible further colli- 
sions in this time step. We feel that these "secondary" atoms would 
not have enough time to have a collision in this time step, and if they 
did, the collisions will be calculated in the next time step. At this point, 
a possibility does exist that some reasonable collisions are excluded 
by the process, but extensive calculations have not revealed any 
examples of this exclusion to date. 

The absolute time of an atom's last collision is also calculated 
during the collision processes described above. For each of the struck 
atoms this time is the total time from the start of the computation to the 
start of this time step plus the time, T, of collision (T is measured from 
the start of the current time step) . 

The time of M's last collision is taken as the start of its collision 
with the last atom on the KHIT list (computed as above) plus a certain 
percentage (PERCENT , an input parameter) of the present time step 
length. M is thus prevented from participating in further collisions 
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for a definite period of time after the start of its last collision. We 
feel that this approach is a valid one , however study of this aspect of 
the model should be continued to determine its true usefulness . 

Sites must now be found for those atoms on the KHIT list whose 
pre-collision energies were above ETHRESH and whose energies are 
now below ETHRESH. 

If the rounded-off co-ordinates of an atom in this class (call it J) 
give the co-ordinates of an actual lattice site , then the atom will 
occupy that site, possibly as part of an interstitial pair if the site is 
already occupied. Atom J will then be allowed to vibrate around its 
site until its energy falls below the THERMAL cutoff. If J's rounded-off 
co-ordinates are not the co-ordinates of a site, then the nearest neighbor 
sites are found and the site closest, space-wise, to atom J is chosen as 
J's vibrational center. Again, this may result in the formation of an 
interstitial pair. Since one of the atoms has an appreciable energy, 
however, this interstitial pair should be destroyed in the next time step. 

The site decided upon may already be occupied by an interstitial 
pair: if so, a "triple interstitial" is formed. An error counter is then 
increased (MISTAKE (5)). The possibility of "triple interstitials" could 
be eliminated by more extensive programming, but an evaluation should 
first be made of the seriousness of this error. The frequency of forma- 
tion of these "triple interstitials" has been small thus far, not more 
than two or three per run, which is tolerable. If this frequency reaches 
an intolerable level in future development, the model must be modified 
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to eliminate any possibility of their occurrence. 



All atoms on the KHIT list whose post-collision and scaling energies 
are below THERMAL are placed on the LATER list at this time. Use of 
the LATER list is discussed below. 

If one of the struck atoms on the KHIT list was a member of an 
interstitial pair, the other atom must now be considered. If this other 
atom was not hit, and the struck atom now has an energy above ETHRESH 
(a very unlikely occurrence) , then the other atom is placed back on the 
site, rather than being left in the split interstitial position. We feel 
that moving the atom back onto the site is more realistic than the 
definite error of leaving it in the split interstitial position. 

All calculations in the model are made on the assumption that the 
atoms are both space-wise and time-wise at the start of the current time 
step, except those that have completed the current time step. Therefore, 
all the atoms on the KHIT list are now moved back in time and space on 
their post-collision tracks to the start of this time step . 

The entire collision process, including the construction of a new 
NUM list is now repeated for atom M. This repetition will continue for 
M until M is unable, time-wise, to have further collisions in this time 
step. The next atom on the LAST list will then become atom M and the 
collision process will be repeated until this M is unable to have further 
collisions in this time step. This overall procedure is then repeated 
until the LAST list has been exhausted. 
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All the atoms on the LAST list, except those on the LATER list, 

(the original M is the first atom on the LAST list) are now moved, space 
and time-wise, to the end of the time step. These atoms have now com- 
pleted the time step and will not enter into the subsequent calculations 
until the next time step . 

The LATER list is now utilized to place the atoms that have partici- 
pated in this sequence of collisions and whose energies are below 
THERMAL on sites, either as members of interstitial pairs or by them- 
selves on previously vacant sites. 

The co-ordinates of the atoms on the LAST list are now examined . 

If an atom has left the lattice, the time step number and the particular 
face (see APPENDIX I, NFACE) of exit is recorded in the LB array. 

This completes the consideration of atom M (the original M) for this 
time step. The next atom (M + 1) is then considered in the same manner 
until all atoms have been through the procedure (and, consequently, the 
time step) . 

After all the atoms have been through these processes , all un- 
affected atoms are moved to the end of the time step, which completes 
the time step. A new time step is then started, or the run is shut down 
by one of the three processes mentioned earlier. 

A numerical example of the order in which atoms are considered may 
be of some value here. Assume atom #73 (M) hits four atoms simulta- 
neously, atoms 75, 79, 87, and 65. Assume further that #73 hits them 
in the order: #79, #65, #87, #75. Atom #73 will hit these four, then 
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repeat the process for further collisions. Assume #73 hits #80 and #67 
next, but has no further collisions after that. Then, atom #79 will 
become the next M since it is second on the LAST list, #65 the next M, 
and after #75 has completed its possible collisions, if any, the regular 
sequence will be taken up once more. Atom #73 was the last one con- 
sidered sequentially, therefore atom #74 now becomes the primary atom, 
atom M . 

As mentioned above , atoms whose energies fall below THERMAL are 
either placed alone on vacant sites , or formed into split interstitial pairs 
with other atoms near them. This is done only after a primary atom, to- 
gether with those on its LAST list, has been completely considered. If 
this were done at the end of every collision process, it would severely 
hamper the multiple collision process, and Impose a great deal of strain 
on the credibility of the argument behind it. 

When a list of the nearest neighbor sites is formed, if the atom in 
question (call it MJ) is within one unit of a face of the crystallite , some 
of the geometrically nearest neighbors normally found will not be present. 
There are 22 special cases for each one of the two general cases (MJ's 
co-ordinates are a site or non-site) . Sixteen of these special cases 
for each concern the edges and comers of the crystal. In our model, 
the 32 (total) possible cases for atom MJ near an edge or comer are not 
considered separately, but are Included in the 12 special cases of when 
MJ is near a lattice face. When one of the 32 unconsidered special 
cases does arise, an error counter is increased. Inspection of all 
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preliminary results to date have, however, shown NO discrepancies as 
a result of this simplification at the edges and corners of the lattice. 

After the list of nearest neighbor sites , the MOON list, is formed, 
these sites are inspected for possible vacancies. If only one vacancy 
exists, atom MJ is placed on that site and the appropriate corrections 
are made to the LB array. If more than one vacancy exists among MJ's 
nearest neighbor positions, then one of these must be selected for MJ's 
occupancy. We feel that this selection should be made on a space-wise 
condition of proximity, rather than on any time-wise condition. There- 
fore, atom MJ is allowed to occupy the closest vacant site. If there is 
more than one vacancy at this smallest distance, MJ will be placed on 
its own site, if available, or the selection will be made in a random 
manner. The choice of a particular vacancy in the case of multiple 
vacancies is somewhat unsophisticated, however at this stage of devel- 
opment in the radiation damage field, we feel that no process of a more 
sophisticated nature can be justified. 

If no vacancies exist among atom MJ's nearest neighbor sites, then 
a prospective interstitial site and partner are chosen. If atom MJ has 
no energy or velocity, this decision is made on a space-wise basis, 
i.e. the closest site is chosen. If atom MJ does have energy, then the 
selection is made on a time-wise basis, i.e. the site that atom MJ 
passes first. 

Once the selection of a prospective interstitial partner has been 
completed, the atom (call it MM) which occupies the site is found. 
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Atom MM's nearest neighbor sites are then found, and the list is 
examined for vacancies. If any vacancies exist, then atom MM is 
moved to one of these vacant sites by the same method discussed 
above for MJ. Atom MJ is then moved onto atom MM’s recently vacated 
site and no interstitial is formed. 

If, however, atom MM finds no vacancies , then an interstitial pair 

must be formed. We feel that of the three types of interstitial pairs 

mentioned earlier, the split interstitial is the most plausible. Indeed, 

Gibson, et al. have shown the split interstitial to be the only stable 
1 7 

configuration. The AX, AY, and £Z distances between atom MJ 
and the site atom MM occupies are found. The largest of these is 
taken as the axis of the interstitial pair. Each atom is then placed on 
its split position, 0.5 units from the site along the axis. Once the 
appropriate corrections are made to the LB array, the formation of the 
pair is completed. 
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7 . Operating Procedures 



The model consists of three separate but inter-connected programs, 
MASTER, SLAVE, and RON. Machine memory limitations forced the devel- 
opment of three , rather than one , programs . 

The main program, MASTER, is the computational program and does 
all of the calculations necessary to complete the run. The two second- 
ary programs, SLAVE and RON, merely use the binary output of program 
MASTER as input and put the results in a semi-analyzed readily legible 
form. The orbits (or intermediate positions of the moving atoms) are 
given by program RON while program SLAVE tabulates and categorizes 
the final positions, energies, etc. of the lattice atoms. Again, be- 
cause of memory limitations, the two secondary programs may not be 
combined . 

If a binary version of program MASTER is obtained with the FORTBIN 
compiler, the storage required may be found by listing of the binary tape. 
The maximum lattice may then be computed . The FORTMAP compiler 
should be used to compile the large version of MASTER since it requires 
less memory than the other compilers. 

Once program MASTER is compiled, a run may be started. MASTER is 
called into memory and the run begins. MASTER calls for BCD input on 
Tape 2, and this is normally transferred to the card reader. BCD output 
is given on Tapes 3 & 4. Tape 4 is normally transferred to the typewriter 
for output (see BOX 10 , APPENDIX V) . The BCD output on Tape 3 is 
temporary and may be eliminated when the program is working correctly. 
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This will also make more memory space available . 

Program MASTER writes binary output on Tapes 6,7, and 8. The 
MASTER output on Tape 6 is the binary input to program RON . The MASTER 
output on Tape 7 is the regeneration output and is used only by the MASTER 
program. The MASTER output on Tape 8 is the binary input to program 
SLAVE. MASTER will automatically rewind each of the three binary tapes 
when output for a particular tape is completed. 

When MASTER has completed computation program SLAVE should be 
called from our system library. SLAVE calls for BCD input (comment 
cards, etc.) on Tape 2, and this is normally transferred to the card 
reader. SLAVE will then read Tape 8, compute, and write its BCD out- 
put on the designated tape. Program SLAVE, at this time, has PRINT 
statements for all output except that mentioned in BOX 1 of the SLAVE 
part of APPENDIX V, and therefore the output will be on Tape 4 unless 
otherwise directed by the operator. 

When SLAVE has completed its run, program RON should be called 
from the system library. The input and output designations for RON are 
the same as those for SLAVE except that RON reads Tape 6 rather than 
Tape 8 . 

If, for some reason, MASTER does not complete its run, and Tape 6 
is not rewound , program RON must be run with Jump Switch 1 in the set 
position. This will insure that the correct END designation (in this case 
a -1) is placed at the end of the information on Tape 6 and that the tape 
is rewound. 
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8 . Parameter Adjustment 



The two input parameters, ETHRESH and TURN primarily govern the 
cohesiveness or solidity of the lattice. At ETHRESH = 3 ev. , for 
instance, some 400-500 replacements occur with a 100 ev. shot. With 
ETHRESH raised to 15-25 ev. the number of replacements drops to con- 
siderably less than ten, depending on the location and direction of the 
initial knock-on atom . The best working values seem to lie in the 
15-25 ev. range, as preliminary results in this range show good correla- 
tion with the work of others . More accurate values of ETHRESH and 
TURN should be obtained in the future from concurrent N-body calcula- 
tions for use in the model. 

Adjustment of the TURN parameter also governs the length of energy 
chains. If TURN =0.0, the energy chains are almost non-existent. 
Atoms vibrating on their sites are then re-directed toward their sites 
after their first collision. Since this first collision is usually a slight, 
glancing one , very little energy is transferred and the energy chain is 
not seen. If, however, the TURN parameter is set too high, displace- 
ments occur among atoms with energy < ETHRESH. Investigations thus 
far seem to indicate that a value of about 0.1 is the optimum value for 
TURN. 

THERMAL has little effect on the larger aspects of radiation damage 
i.e. chains, channels, displacements, ring displacement, etc., but 
has a large effect on machine running time. By changing the value of 
THERMAL from 0.1 ev. to 0.5 ev. , the machine time can be reduced by 
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75% or more . The very low (below THERMAL) energy vibrations are no 
longer seen, to be sure, but since this is not one of the areas of pri- 
mary interest, the advantage of drastically decreased machine time far 
outweighs the inherent disadvantages „ 

When the penetration of higher energy (above a few hundred ev.) 
initial knock-ons is the primary effect under consideration, THERMAL 
may be set as high as 3-5 ev. For initial knock-on energies of 100 ev 
and below, THERMAL should be 0.1 - 0.5 ev. to achieve a reasonable 
degree of thermalization. 

Investigations have shown that the radius of effect (sphere of in- 
fluence) , SPHI, should lie within the range 2. 5-2. 6 Angstroms (for 
copper). Values much outside this range (2.2 say) give rather meaning 
less results. We feel, as a result of the studies we have made, that 

the best value is probably slightly less than the nearest neighbor dis- 

l/2 

tance. For copper, the nearest neighbor distance is (2) 7 (1.807). 

We have used the value (1.414) (1.807) when not specifically investi- 
gating SPHI effects . We feel that any value larger than the nearest 
neighbor distance will produce invalid results. 

Studies of the optimum values for PERCENT and CUTOFF are not 
completed at this time. Preliminary investigations show that values 
above 0.1 for PERCENT and below 0.1 for CUTOFF prevent significant 
collisions from taking place. 

The value of TFAC should be kept below 1.0, unless the normal 
value of 1.0 is used. A more accurate tracing of individual particle 
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movement is possible with smaller TFAC's, but machine time increases 



proportionately. A compromise must then be made between the running 
time available and the individual tracing accuracy desired . 

FMASS may be adjusted for the particular metal being used. At 
present the model is restricted to metals with face centered cubic struc- 
ture . The Born-Mayer potentials used are restricted to copper only 
since values of the interaction input parameters (the constants A & B 
mentioned above) are not available for other materials . 
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9. Preliminary Results 



Preliminary results to date show that the existence of channels, 
energy and displacement chains, vacancies, interstitials, and replace- 
ments depend quite heavily, and quite critically, on the values of the 
input parameters used. Before very satisfactory results can be obtained, 

more reliable values of these parameters must be determined. 

1 7 

A run corresponding to Gibson, et al's. Run No. 12 has been made 
with a slightly modified version of the program, using a 40 ev. primary 
knock-on. The primary knock-on was shot directly at one of the atoms 
on the front face (Gibson et al. start their primary knock-on inside the 
crystallite), at an angle of 15° away from the normal (100) direction. 

The values of the various parameters for the run were: 

1) IX x IY x IZ = 11 x 15 x 8 

2) ETHRESH = 20.0 ev. 

3) THERMAL = 0.5 ev. 

4) CUTOFF = PERCENT = .15 (i.e. 15%) 

5) SPHI =2.5 

6) TFAC =1.0 

7) TURN = 0.1 

Our results (see Figure 1) are similar to those obtained by Gibson, 
et al. in that the energy chains play a dominant role in the dissipation 
of energy. The (100) chain of our run is not as prominent as Gibson, 
et al's. , but strong focusing is observed in the (110) directions. Our 
(100) chain was terminated prematurely by the velocity turn-around 
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Figure 1 . Results of Preliminary Run 
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mechanism employed for lattice cohesiveness . 
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10. Recommendations 



Many ideas and plans occurred to us during the development of the 
model to its present level that we have not been able to institute be- 
cause of time limitations . These unrealized plans and ideas are offered 
as possible guidelines for future development. 

We are convinced at this time that the programming errors in the 
model are concerned almost wholly with the LB information array. As we 
did not anticipate the great dependence upon the LB array when we began 
the problem, we feel that the LB aspects of the entire model (all three 
programs) must be overhauled, and more orderly storage of the informa- 
tion should be investigated . 

Once this overhaul and possible re-arrangement of the LB array has 
been accomplished and the model is working as designed, a more com- 
plete study should be made of the effects of parameter adjustment. 

At present when an atom on the LATER list doesn't find a vacancy, 
a prospective interstitial partner is chosen from the MOON list. A 
search is then made for vacancies near this partner. If none are found, 
the interstitial pair is formed. We strongly recommend that when the 
prospective partner finds no vacancies, the rest of the atoms occupying 
nearest neighbor sites should also be inspected for vacancies near them. 
If, again,, no vacancies are found then the interstitial pair should be 
formed between atom MJ and its original prospective partner. 

To facilitate this , we recommend that the INVAC subroutine be 
split into two or more subroutines . There should be a separate 
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subroutine whose only function is to find the nearest neighbor sites, 
not necessarily by the present system. 

In the present model, if a vacancy occurs near an interstitial pair 
after the interstitial pair has been formed , this Frenkel pair will not re- 
combine unless the interstitial pair is hit by another atom. Provisions 
should therefore be made for periodic inspections of all interstitial 
pairs and any possible vacancies near them. A procedure very similar 
to that discussed above for making interstitial pairs could be used. 

We also feel that subroutines and functions should be incorporated 
to a much greater extent than at present, and that wider use should be 
made of comment cards in the programs. 

A study of the relative merit of a space -wise versus a time-wise 
selection process for the KHIT list has been suggested above and should 
be made. 

The ability to use atoms of two different masses should be built 
into the model. The actual use of two different masses with a Born- 
Mayer potential must be delayed until appropriate constants for the 
potential are found. 

ETHRESH should be made a directional dependent variable rather 
than a constant as at present. Since a function able to represent this 
directional dependence of the displacement energy is not available, a 
table of thresholds associated with direction seems the most promising 
course of action. 
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Two memory saving Ideas have occurred to us: 

1) The possible elimination of the TLAST concept, and 

2) The possible use of a method, already developed, which 
packs the co-ordinates and velocity components of an 
atom into three memory cells rather than the six required 
at present. This sub-program, SQUEEZE, developed by 
one of the authors (RLK) , saves space at the expense of 
running time . 

The use of these two ideas would save considerable memory space, and 
enable the users to approximately double the size of the available 
lattice . 

The energy and velocity scaling processes used after the interaction 
(see BOX 20, APPENDICES III, IV, & V) are a first approximation to the 
problem and should be considered in that light. We feel that at present 
this is one of the weakest areas of the model. Considerable effort should 
be devoted to this section in the future. 
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11. Conclusions 



The model is not fully developed at this time, but we feel that 
some valid conclusions may be drawn from the results obtained thus far. 
The MASTER program (see APPENDIX III), as it exists now, still contains 
some minor programming errors but we feel that the logic is basically 
sound. 

The outstanding attribute of the model at this time seems to be the 
greatly decreased machine time necessary to achieve results compared 
to the machine time necessary for an N -body model. This decreased 
machine time makes it possible to use a lattice as large as the machine 
memory can accommodate . 

We do not feel that this model is sufficiently sophisticated at this 
time to state definitely the degree of confidence that should be placed 
in it. We do feel, however, that the preliminary results show that with 
more work the model can be made into a reasonable approximation to the 
radiation damage problem. We anticipate that our specific recommenda- 
tions, if carried out, will lead to a program which approaches this 
required degree of sophistication. 
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APPENDIX I 



SPECIFICATION OF UNITS AND ANNOTATED LIST OF VARIABLES 
The units used in the model for energy, mass, length, time, and 
velocity are given below: 



ENERGY 


The electron volt (ev.), equal to 1.60207 x 10 ^joules. 


MASS 


The Atomic Mass Unit (AMU) , equal to 1 . 65983 x 10“^ 
kilograms . 


LENGTH 


The Angstrom unit, equal to 10*^ meters. 


TIME 


The jiffy, defined as 10~^ seconds. 


VELOCITY 


Derived. One Angstrom/jiffy = 10^ meters/second. 



The variables given below are essentially in the order of their 
appearance in the program listing. However, certain basic variables 
are given at the beginning of the list. An asterisk (*) denotes an input 
variable . 



VARIABLE 


USAGE 


B (I ,L) 


An array containing the position co-ordinates , velocity 
components, and energies of the atoms in the lattice. 

L is the number of the atom . For I there are seven cases 

1 The X co-ordinate in Angstroms 

2 The Y co-ordinate in Angstroms 

3 The Z co-ordinate in Angstroms 

4 The X velocity component in Angstroms/jiffy 

5 The Y velocity component in Angstroms/jiffy 

6 The Z velocity component in Angstroms/jiffy 

7 The energy in electron volts 
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VARIABLE 



USAGE 



LB (L) 



An array used for information storage , supplementary to 
the B array. L is the number of the atom. Stored from 
right to left we find the following information: 



OCTAL 

POSITION 

1 Bit 1 



Bit 2 
Bit 3 



2-5 



6-9 



10 - 13 



14 



MEANING 

1 if the atom has completed the current time 
step. 

0 if the atom has not completed the current 
time step. 

1 if the atom has an energy > THERMAL. 

0 if the atom has an energy <. THERMAL. 

1 if the site L is occupied by any atom. 

0 if the site L is unoccupied. 

The number of the time step in which atom L 
exited the lattice. Zero if L is in the lattice. 
The number of the other atom in an interstitial 
pair with atom L. Zero if L is not a member 
of an interstitial pair. 

The site occupied by this atom, even if it is 
the original site. If this number is zero, then 
the atom has an energy > ETHRESH and is 
wandering through the lattice. 

The number of the face through which atom L 
exited the lattice, called NFACE. 
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VARIABLE 



OCTAL (L) 
TLAST (L) 



CUTOFF* 



PERCENT* 



USAGE 

OCTAL 

POSITION MEANING 

15 Bit 1 1 if the atom is on the LATER list but has not 
been through INVAC . 

0 if the atom is not on the LATER list. 

Bit 2 1 if the atom has been through INVAC and has 
not had another collision. 

0 if the atom has never been through INVAC , 
or has had a collision since going through 
INVAC. 



Bit 3 Not used 

16 Not used. If Bit 3 is used this results in a 

negative number and should be avoided. 
Equivalent to LB (L) . Used for output purposes only. 

The absolute time from the start of the run to the time of 
atom L's last collision. Given in jiffys. 

A percentage of a Time Step Length. All atoms with a 
time of collision within CUTOFF % of a TSL of the mini- 
mum time will take part in a near simultaneous collision 
with atom L. 

A percentage of a Time Step Length. The absolute time 
of completion of an interaction equals (Time at start) 
plus (PERCENT) (TSL). 
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VARIABLE 



USAGE 



TFAC* 


A scaling factor for the Time Step Lengths . 


LULL* 


The interval, in time steps, at which atoms will be shot 
into the lattice, provided the proper jump switches are 
set. Does not affect operation if the jump switches are 
not set. 


IX, IY, IZ* 


The X, Y, and Z fixed point dimensions of the lattice, 
respectively, such as 9 x 14 x 15. IX must be ODD. 


A* 


One half the side of a unit cell. For copper this is 
1,807 Angstroms . 


SPHI* 


The sphere of influence of an atom. One atom has an 
SPHI of (1.4 14) (A), and the other atom is considered 
to be a point. Given in Angstroms. 


ETHRESH* 


The displacement threshold energy in electron volts. 
An atom with an energy > ETHRESH may move through 
the lattice. Atoms with energy below ETHRESH are 
forced to remain on their sites in the lattice. 


THERMAL* 


The thermal energy limit. Atoms with energy above 
THERMAL are allowed to vibrate on their site. Atoms 
with energy below THERMAL are placed on site? , or in 
interstitial pairs by the INVAC subroutine . 


FMASS* 


The mass of a lattice atom in Atomic Mass Units (AMU) . 


JPB* 


The interval, in time steps, at which the Time Step 
Length is recalculated. 
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VARIABLE 


USAGE 


MNTS* 


The maximum number of time steps the program will run. 


MNP* 


The maximum number of particles to be shot into the 
lattice . 


FI (1-3) 


Used throughout the program. Usually denotes the 
floating point co-ordinates of a particular site. 


AB (1-3) 


The physical dimensions of the lattice in Angstroms in 
the X, Y, and Z directions, respectively. 


NLINE 


The number of atoms in a line in the X direction for the 
undisturbed lattice . 


NPLANE 


The number of atoms in any undisturbed X - Y plane of 
the lattice. 


LMAX 


The maximum , or total , number of atoms . 


MAXLAT 


The total number of atoms in the lattice , excluding 
those shot into the lattice, i.e. the number of sites. 


TEMP, TEM 


Floating point variables used for temporary storage only, 


SITE (L) 


A short subroutine. Given the number of a site as an 
input, SITE computes the fixed point co-ordinates of the 
site. These are then found in IB (1-3). 


TTIME 


The total elapsed time from the start of the run to the 
end of the current time step, in jiffys. 


OENERGY 


The total original energy of the lattice in ev. This is 
increased every time a particle is shot into the lattice . 
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VARIABLE 



USAGE 



ENERGY The energy in ev. of the most energetic atom in the 

lattice. Atoms that have left the lattice are not con- 
sidered in the computation of ENERGY. 

TENERGY The total energy in ev. in the lattice. This is computed 

every time step for comparison with OENERGY. 

IENERGY The number of the most energetic atom in the lattice . 

MISTAKE (l — 1 0) An array of 10 error counters. Their meanings are: 

1 The number of time steps for which TENERGY varies 
more than 1% from OENERGY. 

2 The number of times a vacancy is filled, or an 
interstitial is formed, near an edge or corner of 
the lattice. 

3 The number of times a negative energy was computed. 

4 The number of times an atom formed an interstitial 
with itself, i.e. a false interstitial. 

5 The number of triple interstitials formed. 

6-9 Not used . 

10 Set equal to MISTAKE (2) at the start of INVAC so 
that the increase in MISTAKE (2) may be calculated. 

NA The index for the TSLI and NTSC lists . Always one 

greater than the number of times the TSL has been 
recomputed . 
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VARIABLE 


USAGE 


NEXT 


The number of the next particle to be shot into the 
lattice. NEXT is one before any atoms are shot into 
the lattice. 


MINTS 


The minimum number of time steps, i.e. the lower limit 
of the time step DO loop. Always one unless the re- 
generation method is used. 


Til 


The square of SPHI. 


MOVE 


An indicator used to determine the exit point from sub- 
routine INVAC . 


NOW 


The lattice site nearest an atom, computed in INVAC, 
used to determine an atom’s oscillatory center, if the 
atom' 3 energy is below ETHRESH. 


LOCATE 


An indicator used to hold the number in Index Register 5, 
enabling us to use Index Register 5 without destroying 
the contents. 


CON (1-4) 


Constants used in the subroutine INTER. These have no 
relation to the main program. 


W (1-4) 


Constants used in subroutine INTER. These have no 
relation to the main program. 


XC (1-6)* 


Only two used. These are the constant parameters used 

in the Born-Mayer potential , and can be changed to fit 

any one of the three types of Born-Mayer for copper, 

17 

as devised by Gibson, et al. 
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VARIABLE 


USAGE 


Cl - CIO 


Constants used In subroutine INTER. Computed outside 
the subroutine as a time saving measure. 


MASK1 - 
MASK8 

\ 


Masking quantities used for logical arithmetic with 
LB (L) o The definition of the significance of the octal 
positions in LB (L) and the masks themselves explain 


ALFA 


their specific applicability. 

10~® A small number for comparison to floating point 
variables to prevent use of insignificantly small quanti- 
ties . 


ALPHA* 


The angle between the original X and Y velocity compo- 
nents of a particle shot into the lattice, given in radians. 


BETA* 


The angle between the original X and Z velocity compo- 
nents of a particle shot into the lattice , given in radians . 


YENTRY (NEXT) 


The Y co-ordinate at which a particle shot into the lattice 
crosses the X = 0 plane . Given in Angstroms . 


ZENTRY (NEXT) 


The Z co-ordinate at which a particle shot into the lattice 
crosses the X = 0 plane. Given in Angstroms. 


EPB 


The energy in ev. of a particle shot into the lattice. 


TSL 


The Time Step Length in jiffys. TSL is calculated so that 
the most energetic atom in the lattice will move a dis- 
tance equal to (TFAC)(A) in one time step. 


TSLI (NA) 


A list of the first 30 TSL's that are calculated. Used 
for output purposes only. 
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VARIABLE 


USAGE 


NTSC (NA) 


A list, corresponding to the TSLJ list, that gives the 
time steps in which the values given in the TSLI table 
were calculated. 


EMAX 


A list of the values of ENERGY for the first 30 time 
steps. Used for output purposes only. 


N 


The sixth Index Register, and the number of the current 
time step. The current time step can always be found 
by stepping the computer console and recording the 
value of this Index Register. N is used for nothing else. 


L 


The number of the atom sequentially under consideration. 
The index of the DO loop starting at Statement 160. 


M 


The number of the atom currently under primary considera- 
tion, usually equal to L. 


J 


Usually the number of the atom considered in relation to M 


LATER (NZ) 


A list of atoms whose energies are less than THERMAL 
and are to be put through the INVAC subroutine . 


NZ 


The index for the LATER list, used for no other purpose. 


INMIN 


The index for the LAST list, used for no other purpose. 


INMAX 


The maximum number of atoms listed on the LAST list, 
equal to NMAX when the list is made. 


INPUT 


An indicator. If zero, then the KHIT list is copied onto 
the LAST list. If not zero, the LAST list is left untouched. 
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VARIABLE 


USAGE 


NUM ( ) 


A list of all the atoms in a box of side length 3.02 (A) 
centered at atom M . Only atoms on this list are con- 
sidered as possible targets for M. M is first on the 
list. 


KMAX 


The total number of atoms listed on the NUM list. 


KHIT ( ) 


A list of atoms, derived from the NUM list, that atom M 
will hit. Atom M is first on the list. 


NMAX 


The total number of atoms listed on the KHIT list. 


TMIN ( ) 


The time, in jiffys, of closest approach between atom M 
and the corresponding atom on the NUM list. This 
assumes no deviation from straight paths. Measured 
with the start of the current time step as a zero reference 


DSTANCE ( ) 


The distance, in Angstroms, between atorrts when at their 
point of closest approach as determined by their TMIN. 


T ( ) 


The time, in Jiffys, at which the two atoms are SPHI 
distance apart. Will be less than the corresponding 
TMIN if there is to be a collision between M and the 
atom on the corresponding NUM list. 



ITEMP, JTEMP , Fixed point variables used throughout the program for 
MTEMP, NTEMP 

temporary storage purposes. No particular significance. 



NPAIR 


An indicator for a computed GO TO statement. NPAIR is 
3 if one of the atoms being hit is a member of an inter- 
stitial pair. 



58 



VARIABLE 


USAGE 


IN3 


An indicator that shows the number of those on the NUM 
list that will be hit by atom M . 


LAST ( ) 


A duplicate of the first KHIT list comprised for M when 
M = L. All atoms on the list will become M's when the 
preceeding M is no longer able to have collisions in this 
time step. 


T1 - T10 


Temporary variables (see BOX 17, APPENDIX IV) used to 
speed up computation time . 


NFACE 


The face of exit of an atom when it leaves the lattice. 
Faces 1,3,5 are the negative X, Y, and Z faces, res- 
pectively, while faces 2,4,6 are the positive X, Y, 
and Z faces respectively. 


JA 


The index for the KHIT list that shows which atom 
atom M is colliding with at the present time. JA's range 
is from 2 to NMAX. Has other uses in INVAC . 


KMIN 


The number on the T list that is smallest. T (KMIN) is 
the minimum time, and is the collision time for atom 
NUM (KMIN). 


NHIT 


The number of atoms on the KHIT list that have been 


ELOST 


hit by M. 

The energy in ev. that M lost in its collisions with the 
other atoms on the KHIT list. 
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VARIABLE 



USAGE 



SAVE (l,K) 



COLD ( ) 



PAR 



TIME 



DEV (1,1-2) 



The pre-collision velocity components and energies of 
the atoms on the KHIT list. SAVE (1,1) is the initial 
X velocity of KHIT (1) . SAVE (2,1) is the initial Y 
velocity of KHIT (1) , etc. SAVE (4,1) is the energy of 
KHIT (1). 

A floating point variable used in the energy and velocity 
scaling process. COLD (1-6) correspond to the first 
six co-ordinates in the B array. They are the sums of 
the individual components of the DEV array after each 
collision, i.e. COLD (1) is the sum of all the DEV 
(1,1) 's for a series of collisions . COLD (1) is then 
the sum of all the X co-ordinates M would have had as 
a result of the collisions with the rest of the KHIT list. 
The scattering parameters in meters , and the third input 
to the INTER subroutine. 

The time T for M's previous collision with one of the 
others on the KHIT list. For instance, if JA is 3, we 
are considering KHIT (3) in relation to M„ TIME is then 
the T corresponding to KHIT (2) . 

An array, used in subroutine INTER, holding the position 
co-ordinates, velocities, and energies of the two atoms 
at the end of the interaction. The first index, I, has 
the same meaning as the first index of the B array. 
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VARIABLE 



USAGE 





The second Index is one for atom M, and two for atom J. 
The B array is not changed in the subroutine , but the 
DEV array contains the output of the subroutine. The 
B array is then scaled as necessary. 


STUPID 


A scale factor used in the velocity scaling process. 


SCALE ( ) 


A three member array used in the velocity scaling process. 


TURN* 


An atom in the intermediate energy range must lose TURN % 
of its energy before it is re-directed back toward its 
oscillatory center. 


IT 


A temporary variable used to facilitate programming. Also 
the index for the DO loop over the LATER list. 


MANY, NANY 


The atoms occupying the site chosen for atom M after M's 
energy has fallen below ETHRESH. If NANY is zero, then 
M and MANY will become an interstitial pair. If both are 
non-zero, a "triple interstitial" is formed, while if both 
are zero, M occupies the site alone. 


JPMAX 


The number of atoms on the MOON list, i.e. the number of 
nearest neighbor sites. 


MOON ( ) 


A list of nearest neighbor sites for a particular atom. 


NOW1 , NOW 


Temporary variables used to indicate the site chosen as 
oscillation center for an atom in the intermediate energy 
range. NOW1 = (NOW shifted 27 binary positions to 
the left) . 
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VARIABLE 


USAGE 


NAVY, MAD 


Similar to LOCATE, except these are used for Index 
Registers 2 & 1, respectively. 


ZZ (6) 


The time in seconds necessary for the interaction. 


NN, IL, LL 


Variables used to denote atom numbers. NN is the 
index for the DO loop starting at Statement 241. IL 
usually denotes a site number, while LL usually de- 
notes the number of the "other" atom in an interstitial 
pair . 


MJ 


The number of the atom used as an input to INVAC . 


MM 


Usually equals MJ, but if MJ finds no vacancies then 
MM is the prospective interstitial partner atom. 


I FACE ( ) 


A three member array denoting the face an atom entering 
INVAC is near (i . e . within one lattice unit) . The index 
is 1, 2, or 3 for the X, Y, or Z directions, respectively. 


NOON ( ) 


The list of vacancies near atom MM (which is usually MJ) 
This list is derived by examination of the MOON list, 
i.e. if an atom is not on the MOON list, it will never 
be on the NOON list. 


MAX 


Either 1 , 2, or 3. Denotes the direction (X, Y, or Z re- 
spectively) of the axis of the prospective interstitial 
pair. 



NOTE: This is the end of the list of variables for the MASTER program. 
Many of the temporary indices are not listed in this table. 
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The symbols listed below are peculiar to program SLAVE and the 



usage given applies only to program SLAVE: 



VARIABLE 


USAGE 


NTS 


The actual number of time steps the run completed. 


KET (1-26) 


A list of the energy ranges for the histograph distribu- 



tions . 



IHISTO (1-105) 


The histograph format array. 


NBUL 


The number of knock-ons (i.e. atoms with energy above 



ETHRESH) in the lattice „ 



EBUL 


The total energy of ail knock-ons in the lattice. 


NSUBST 


The number of replacements in the lattice. 


NSAME 


The number of atoms on their original lattice sites. 


NVAC 


The number of vacant lattice sites . 


INTERST 


The number of interstitial pairs in the lattice. 


ETHERM 


The total energy of the atoms that have gone through 



INVAC. 



NSIDE (1-6) 


The total number of atoms that have left the lattice 



through each of the six faces . 



ESIDE (1-6) 


The total energy of the atoms that have left the lattice 



through each of the six faces. 



ECAL 


The same as TENERGY, except calculated by Program 



SLAVE. 



LMAXCAL 


The same as LMAX, except calculated by Program SLAVE 


YANGLE 


Angle ALPHA given in degrees rather than in radians. 
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VARIABLE 


USAGE 


ZANGLE 


Angle BETA given in degrees rather than in radians. 



NUMBERS (1-50) A fixed point array used for temporary storage only. 



IC (1-100) 


A fixed point array used for temporary storage of the 
histograph channels . 


DEPTH 


The total penetration of an atom in the positive X 
direction . 


PERDIST 


The total radial penetration of an atom perpendicular 
to the X direction . 


LL 


Usually the number of an atom, temporary usage only. 


XYZ (1-3) 


Corresponds to the FI array in program MASTER. 


NTOTAL 


The sum of all the numbers on the NUMBERS list. 


NUMBIG 


The largest number on the NUMBERS list. 


TOTALN 


NTOTAL in floating point format. 


BIGNUM 


NUMBIG in floating point format. 



The symbols listed below are peculiar to program RON and the usage 
given applies only to program RON: 



VARIABLE 


USAGE 


T 


Temporary floating point variable. Not to be confused 
with the T array in program MASTER. 


DX 


A distance parameter, in Angstroms. If any of an atom's 
co-ordinates differ from its stored values by DX, the 
atom's position will be printed, and the new position 
will be stored . 



64 



* a' tii 


















VARIABLE 


USAGE 


RA 


Reciprocal of A. 


P (1-3) 


The co-ordinates of the atom under consideration. 


IN 


The number of the time step the atom under considera- 
tion has co-ordinates P (1-3). 


X (1-3 ,N) 


Stored co-ordinates of an atom to be written out at a 
later time. 


NT (I) 


An array used in the output section. The lower half of 
a storage cell is the number of the atom whose co- 
ordinates are X (1-3, N) , and the upper half is the time 
step the atom had these co-ordinates. 


NO 


An indicator. If greater than zero, a heading has been 
written on the output tape for the atom under considera- 
tion. If less than zero, no heading has been written 
on the output tape for the atom under consideration . 


S (1-3) 


Co-ordinates of an atom in floating point form in A units . 
One A unit is 1.807 Angstroms for copper. (See A, above) 
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APPENDIX II 
LATTICE CONSTRUCTION 




0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 



FIG. II (100) LATTICE PLANE, Z EVEN 
>Y IX = 15 IY = 14 
NLINE = (15+1) /2 = 8 

NPLANE =(NLINE) (IY+1) = (8) (14+1) = 120 
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x 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 

T FIG. Ill (100) LATTICE PLANE, Z ODD 

I y 

To find the original co-ordinates of atom 77, say: 

Z =(77-1) /NPLANE = 76/120 = 0 Remainder = 76 

Y =(76) /NLINE = 76/8 = 9 Remainder = 4 

X =(4) (2) = 8 

Y + Z = 9 + 0 which is ODD, therefore X = X + 1 = 9 
The original co-ordinates of atom 77 were then (7,9,0). 



67 




' 

- 

rX ' • L u • 9 ^ 

I 1 If * m> i jB 5 

■ 



















. 






• . 









. ■ 







■ 



, 






APPENDIX III 
PROGRAM LISTINGS 



PROGRAM MASTER 

DIMENSION 8(7 ,1500) ,IB( 15CC) .OCTAL (1 SOC) , FI ( 3) . IP (3 ), T( 20) 
DIMENSION TMIN120) , IF ACE (3) .MISTAKE! 1 0 ) , MOON ( 1 5 ) ,NOON( 10),W(4) 

B1SIKII8S ,c, 

DIMENSION ZENTRY( 10) , TSL I ( 30 > , NTSC ( 30 ) , TL AS T ( 1SCC) ,LAST(2C ) 
DIMENSION EMAX13C) ,COLD( 7) ,SAVc(4. 10 ) , SCALE! 3) 

COMMON 6, LB, OCTAL, FI , IB,T, TMIN, I FACE , M I STAKE , MOCN, NOON ,W , XC , XP , Z Z 
COMMON CON.OEV, AB, A,NPLANE,NLINE ,LMAX.MASK1 , M ASK2 , M ASK3 , M A SK4 
COMMON I TEMP, JTEMP , MTEMP ,NTEMP , I N2,I N3, IL ,JA , JB, JC, JK, JL, JM, JP, JR 
COMMON JT, JW, MM, KM IN, MAX ,TEMP,NI , JPMAX,NIF, I X. IY. I Z , NN.M A XL A T , JJ 
COMMON IP,TE,TEM,C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C1 1,C12, T1 ,T2 
COMMON DSTANCE,FMASS, 13, TSL 
COMMON MASKS. MASK6,MASK7 , MASKS, MOVE, NCI* 

EQUIVALENCE ( XP » TM I N ) , (LB, OCTAL) 

BOX 1 

WRITE OUTPUT TAPE 3, 995C 
9950 F0RMAK2UH THIS IS PROGRAM MASTER /) 

DO 9980 1=1,5 

READ INPUT TAPE 2,9981 

9980 WRITE OUTPUT TAPE 3,9981 

9981 FORMAT ( 80H 

1 ) 

BOX 2 

READ INPUT TAPE 2 , 803 , CUTOFF , PERCENT , TF AC , TURN 
803 FORMAT! 4F 1 0, 6 ) 

RE AO INPUT TAPE 2 , 1 00 1 , LLLL , I X , I Y, IZ , A 
1001 F0RMAT(4I4,F10.6) 

READ INPUT TAPE 2 . 1 02 , SPH I , E THRE SH ,T HERM AL , FM ASS , JPR, MNT S , MNP 
102 F0RMAT(4F10.6,3I5) 

BOX 3 



FI ( 1) = IX 
F I ( 2 ) = IY 
F I ( 3 ) = IZ 
DO 1000 1*1,3 
1000 AB( I ) =A*F I ( I ) 



BOX 4 

LOA( IX) , INA( 1 ) ,ENQ(0) , ARS( 1 ) ,STA(NLI NE ),LDA( IY) . INA( 1) , 
STA ( NPLANE ) , LOA ( I Z ) , I NA ( 1 ) ,-MUI (NPLANE) ,STA( LMAX) ,STA( MA 



, MU I ( N L I N E ) 
STA( MAXLAT ) 



BOX 5 

DO 1020 L= 1 • LMAX 
CALL SITE(L) 

1010 DO 1 C 1 5 1*1,3 
TEMP* IB ( I ) 

1015 B( I , L ) *A* TEMP 

ENA4(0B) ,ALS(27) ,A0D(48) , ST A4 ( LB ) , LD A ( - 1 0. 0 ) , ST A4 ( TL AST ) 

00 1C20 1*4,7 

1020 B ( I , L ) = 0.0 

BOX 6 

WRITE TAPE 6, I X , I Y , I Z , AD . MAXL AT . A , SP F I , FMAS S , ETHRESH , THE RM AL , MNP 
WRITE TAPE 6 , ( ( R ( I , J ) 1 = 1,3) J* 1 , M AX L AT ) , NL I NE , NPL ANE 
LOA(C.O) , STA (TTI ME), STA(CENERGY) , STA (ENERGY) , ST A ( Z Z* 4 ) , EN A ( C ) 

S T A ( M I STAKE* 1 ) , STA ( M I STAKE + 2 ) , ST A ( MI STAKE + 3 ) . ENA ( 1 ) ,STA( NEXT) 
STA(NA), STA (MINTS) ,LDA( SPHI ) ,FMU(SPH I ) ,STA( T 1 1 ) 

ENA ( 0 ) • STA ( M I STAKE + 4 ) , STA ( M I ST AKE + 5) , ST A ( MOVE ) , ST A ( NOW ) 

STA ( LOCATE ) 
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non oooooonoo 



ecx 7 



START SET UP FOR INTERACTION SUBROUTINE 
CONI 1 )=4. 820781099625 E-02 
C0N(2)=1 .08906256893SF-01 
CON ( 3 1=4.4888 72 968795E-0 1 
CONI 4 1*8. 6595709 11259E-01 






CC 



50 



LDAI W + 2) ,STA(W+3),LDA(W+1 ),STA(W+4) 

READ IN POTENTIAL PARAMETERS 
READ INPUT TAPE 2, 50 , ( XC ( I ) 1=1,6) 

F0RMAT(10X,6F10.5/) 

Cl=l .0E-10*SPHI 
C2=XC( 1 )* 1.6021 E-l9*EXPF (XC(2) 1 
C3 = XC 1 2 1 /C 1 

C4=l .6598E-27«FKASS*0.25 
C5= 1 . OE-1 4 /C 1 
C6=0 . 5E 1 0*C 1 
C 7=C 1 *C 1 
C8= 1.0/Cl 

C9=XC ( 1 1 *1 .6021 E-19 
Cl 0=1 .0/1 FMASS*0. 51 80076 15CC1- 
C ENO SET UP FOR INTERACTION SUBROUTINE 

BOXES. 8 AND 9 

MASK 1=7777000000 
MASK2=77770000C00008 
MASK2=77777777C0000B 
MA SK 4=7000000 00 77770B 
MASK5=1 

MASK6=100000000000000B 
MASK7=2*MASK6 
MASK8=7677777 7777 77 7770 
ALFA=1.0E-6 

USE OF SENSE SWITCHES SELECTIVE JUMPS) 

UP, UP, UP READ BINARY INPUT OF INTERRUPTED PROGRAM 

UP, UP, DOWN WRITE BINARY CUMP ON TAPE 7 FOR REGENERATION 
UP, DOWN, DOWN READ IN NEXT PARTICLE IF THE TIKE IS RIGHT 
DOWN, UP, DOWN JUMP TO 8001 WRITE OUTPUT AND RETURN TO 116 

DOWN, DOWN, UP AT 116, I F SAT I SF I ED , GO ON TO ICO, OTHERWISE JUMP TO E ND , E N 
UP. DOWN, UP WILL WRITE TAPE 8 8 1 NARY, BUT NOT THE BCD PART OF OUTPUT. 
WHILE PRINTING OUTPUT, IF WE PUT STOP SWITCH 1 UP, IT WILL STOP AT THE EN 
OF THE OUTPUT TO ENABLE LS 10 RESET THE PROPER SW I TCHES, ETC . 

BOX 10 

SL J 1 ( 80 1 » SL J C 9 1 1 

80 SLJ2I81 1 , SL J (91) 

81 SLJ3I 90) , SL J ( 91 ) 

90 READ TAPE 7,M I NTS , MNT S , L M AX , B , LB , TLA 
1TSLI ,NTSC,TT I ME, MI STAKE, YENTRY,ZENTR 
2SPHI , FMA SS»JPB»MNP,LAST,tMAX,VMCK,CO 
REWIND 7 

DO 9000 N=MI NTS, MNTS 

WRITE OUTPUT TAPE 3,8 100 ,TENERGY,N,I 
WRITE OUTPUT TAPE 4 , 8 1 00 , TE NERGY , N ,1 
FORMAT! FI 5. 9, 2 1 8, FI 5. 9) 

I ENERGY=0 

THIS IS THE START OF THE MAIN DO LOO 
STEPS. WE CAN CHECK THE CURRENT TIM 
NOTING THE CONTENTS OF INDEX REGISTE 
SL J 1 (82), SLJ (93) 

SL J 2 ( 83 ) , SL J ( 93 ) 

SLJ 3 ( 93 ) 

WRITE TAPE 7,N,MNTS,LMAX,B,LB,TLAST, 

1TSLI , NT SC, TTIME, MISTAKE, VENTRY.ZENTR 
2A,SPHI,FMASS, JPB, MNP.LAST, EMAX, VMCM, 

REWIND 7 
93 CONTINUE 



ST, TENERGY.C ENERGY, E NERGY, TSL, 
Y , NEXT , NA , EPB, LULL, I X, IY, IZ, A, 
LC, EFOUND 



91 



8100 



82 

83 

92 



ENERGY, ENERGY 
ENERGY, ENERGY 



P CF THE PROGRAM, THE ONE ON _ 
E STEP BY STEPPING THE 1604 AND 
R 6. 



TIME 



TENERGY.CENERGY, ENERGY, TSL, 
Y, NEXT, N A, EPB, LULL, IX, IY, IZ, 
COLD, EFOUND 
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BOX 1 1 



SLJ 1 ( 100) ,SLJ2(85) ,SLJ( ICO) 

85 SLJ3( 100) ,SLJ(80C1 ) 

116 SLJ1 ( 9999 ),SLJ2( 9999) » SL J3 ( ICO) , SLJ( 9999) 

100 LDA (NEXT), INA(-1),AJP(12C ),SLJ 1(84), SLJ ( 110) 

]?0 !/iji>(800l ) ,A JP2( 1 30) , SLJ (8001 ) 

115 EN A6 ( 1 8 ) • DV l (LULL),0JP1( 130) , LDA ( NEX T ) , SUB ( MNP ) , A J P( 1 20) ,A JP2 ( 13C ) 

120 RAO(LMAX) ,LIL1 ( LMAX ) , ENA ( 2 ) , ST A 1 ( LB) ,LDA(-2 .0), STAl( TLAST) 

900 REAO INPUT TAPE 2 , 500 , ( 8 ( I , L MAX ) I =1 , 3 ) , B ( 7 , LM AX ) , ALP HA, BET A 
500 FORMAT( 6F 10.6 ) 

YENTRY(NEXT)=B(2,LMAX) 

ZENTRY(NEXT)=B(3,LMAX) 

RAO( NEXT ) 

ZZ( 1 )=TANF( ALPHA ) 

FMU ( ZZ+ 1 ) STA ( ZZ +2 ) 

ZZ( 3 I = T AN F ( BETA) 

FMU ( ZZ+3 ) FAC ( ZZ+2 ) FAC(l.O) STA(ZZ+2) 
B(4,LMAX)=SGRTF(B(7,LMAX)*C10/ZZ (2)) 

B(5,LMAX)=B(4,LMAX)*ZZ( 1 ) 

B(6»LMAX)=B(4,LMAX)*ZZ(3) 

B(3,LMAX)=B( 3,LMAX)-B(6, LMAX ) /B ( 4 , LM AX ) 

B(2,LMAX)=B(2,LMAX)-B(5, LMAX)/R(4,LMAX) 

125 0ENERGY=0ENERGY+B(7»LMAX) 

TENERGY=TENERGY+B(7»LMAX) 

EPB = B ( 7, LMAX ) 

TSL=A*0.71 97274 588/SORTF ( ABSF ( B ( 7 . LM AX ) /FMASS ) ) *TF AC 
+LDA ( NA ) , INAl-30) , +AJP ( L+ 1 ) , A JP2 ( 1 40 ) , L I L 1 ( NA ) , LDA ( TSL ) 

ST A 1 (TSL I ) ,ENA6(0) ,STA1 (NTSC) , RAC (NA ) ,SLJ( 140) 

BOX 12 

130 ENA6( 0) , ENQ( 0 ) ,OVI ( JPB ) , CJP1 ( 1 40 ) 

135 TSL=A«0. 71 97274588/SORTF (ABSF (ENERGY /FMASS) )«TFAC 
+LD A ( NA ) , INA(— 30) ,+AJP(L+l ) ,AJP2( 140 ) 

LI LI (NA) ,LDA( TSL) ,STA1 (TSL I ) , ENA6(0) , ST A 1 ( NTSC ) , RAO ( N A ) 

BOX 13 . 

140 LDA ( TTI ME ) »FAD(TSL)»STA( TTIME) 

LDA ( OENERGY ) ,FMU( 1 .01 ) ,FSB( T ENERGY) , A JP ( 155 ) , AJP3 ( 1 50 ) ,LDA ( 0 .C ) 

FSB (OENERGY) .FMU ( 0.99.) , FAD ( T ENERGY ), A JP ( 1 55 ) , A JP2 ( 1 55 ) 

150 RAO ( M I STAKE+ 1 ) 

BOX 14 

155 EN A6 ( 0 ) , I NA (-3 1 ) , +A JP2 ( L + 3 ) , LO A ( ENER GY ) , I N I 6 ( - 1 ) , STA6 ( EM AX ) 

♦ INI6( 1 ) ,+LDA(0.0) ,STA( ENERGY) ,STA(T ENERGY) 

DO 160 L=1,LMAX 
TENERGY=TENERGY+B(7,L) 
fLDA4(LB) ,SCL( MASK5) , ST A4 ( LB ) 

160 CONTINUE 

BOX 15 

C THIS IS THE START OF THE MAJOR DO LOCP CN THE PARTICLES 
DO 8000 L= 1 , LMAX 

WITHIN A PARTICULAR TIME STEP, THIS IS THE START OF THE MAIN CC LOOP OV 
ALL THE PARTICLES IN THE LATTICE. WE CAN TELL WHICH PARTICLE IS BEING 
SIOEREO BY STEPPING THE 16C4 AND CHECKING INCEX REGISTER 4. 

LDQ4 ( LB ) , LOL ( 2B ) , A JP ( 8C0C ) , L DL ( 77771 E1.AJP1 (8000) 

LDL (MASK7)»AJP1( 8C00) 

165 DO 161 K= 1 » 50 

161 LATER( K )=C 

ENA ( 1 ) ,STA(NZ) ,ENA( 2) ,STA( IN MIN) ,ENA (0) ,STA( INPUT ) 

DO 164 K=1 ,20 

ENA(0),STA3( NUM ) ,STA3(KHIT),STA3(LAST) , LCA(C.O) ,STA3(T) 

STA3 ( DST ANCE ) ,STA3(TMIN) 

164 CONTINUE 

. S I L 4 ( ITEMP),LIL5( ITEMP) 

C THIS IS THE END OF THE USE OF L AS SUCH FOR THE REST OF THE CC LCCP 
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BOX 16 



200 ENA(-C) v STM IN 3) , ENA5 ( 0 ) , A JP ( 277 ) , ENI 3 ( 2 ) ,SIL5(NUM*1 ) ,SIL5(KHIT+1) 
L0Q5 (LB) ,LDL(MASK7) . A J P 1 (277) 

NPA I R= 1 
DO 201 1=1,3 

201 T(I1=B(I,M) 

LD02?LB ) * A JP 1 ( 20 3 ) 

DO 2C2 1=1,3 

IF(ABSF(B(I»J)-TII ) > - 1 . 5 1 * A ) 202,20 3,20 3 

202 CONTINUE 

SIL2( I TEMP) ,LAC( I TEMP) , I NAS ( 0 ) ,AJP(203) ,ENA2(0) ,STA3(NUM ) ,INI3( 1 ) 

203 CONTINUE 

I NI 3 (- 1 ) , ENA3 ( 0 ) ,STA(KMAX) , INA<-2) ,A JP3(277) ,ENA<0) ,STA( IN3) 
STA(NFACE) 

BOX 17 

DO 215 K= 1 , KMAX 
LDA3 ( NUM ) ,STA( I TEMP ) 

T 1 = B ( 1 , I TEMP )-B ( 1 , M ) 

T2=B(2,ITEMP)-B(2,M) 

T3=B( 3, ITEMP)-B(3,M) 

T4=B ( 4 , ITEMP )-B( 4,M) 

T5=B(5, ITEMP )-B( 5 »M ) 

T6=B(6, ITEMP) -B (6, M) 

T7=T4*T4+T5*T5+T6*T6 
T8 = T1*T1 +T2*T 2+T 3* T3 
LDA ( T7) , FSB (ALFA) ,A JP3(2 14 ) 

T9=T1*T4+T2*T5+T3*T6 

LDA ( 1.0) , FOV ( T7 ) , STA ( T 10 ) , L AC ( T9 ) , FM U ( T 1 0 ) , S T A 3 ( TM I N ) 
DSTANCE(K)=SQRTF(ABSF( (T1+T4*TMIN(K) )«(T 1 + T4*TMIN(K) ) + (T2 + T5*TPIN( 
IK) )«(T2+T5*TMIN(K) )+(T3+T6»TMIN(K) ) * ( T 3+ T6* TM IN ( K ) ) ) ) 

F SB ( SPH I ) ,AJP(214) ,AJP2( 214) 

T(K)=TMIN(K)-SQRTF(ABSF(TMIN(K)*TMIN(K)-(T8-T1 1 )/T7) ) 

211 LD A3 ( T ) , FSB ( TSL ) , AJ P( 2 1 4 ) , A JP2 ( 2 14 ) . FAC ( TSL ) ,F AD ( TSL ) , AJ P3 ( 2 1 4 ) 
L0A3(T) ,FDV(TSL) ,FSe( 0.95999 ) ,AJP( 21 2) ,AJP2(214) 

212 LIL1 ( ITEMP) ,L0A3(T) , FAO ( TT I ME ) ,FSB ( T SL ) , FSB 1 ( TLAST ) , AJP( 2 1 4 ) 
AJP3(214) ,LDA3(T) , FAD ( TT IME ) , F SB ( TSL ) , FS B5 ( T L AST ) , A JP ( 2 1 4 ) 

A JP 3 ( 2 1 4 ) 

213 RAO ( IN3 ) , SL J ( 2 1 5 ) 

214 LDA(TSL) ,FMU( 1 . 0 E+2 ) , STA 3 ( T ) 

215 CONTINUE 



BOX 18 

LDA( IN3) »AJP1(216),SLJ(2?7) 

216 LOA (KMAX) .INA(-1 ) , STA( JA ) ,LCA( TSL) ,FMU( 1 .CE+3) , STA (TEMP) 
00 217 K=2 , KM AX 

^ LDA3 ( T ) ,FSB( TEMP) ,AJP2 (217), LDA3 ( T ) , STA ( TEMP ) , SI L3 ( KM IN ) 
2fl7 CONTINUE 

L I L 1 (KMIN) ,ENI2( 2) ,LDA (TSL) ,FMU( CUTO FF ) , STA ( T1 ) 

00 219 K= 1 ,KMAX 

LDA3 ( T ) , FSB 1 ( T ) , F SB ( T 1 ) , A JP ( 2 1 8 ) , A JP 2 ( 2 1 9 ) 

218 L0A3 ( NUM ) ,STA2(KHIT) , INI21 1 ) 

219 CONTINUE 

I N I 2 ( — 1 ) , SIL2(NMAX) ,LDA( INPUT) ,AJP1( 301 ) ,SI L2( INMAX) 

DO 300 K= 1 , INMAX 
LDA3(KHIT) ,STA3( LAST) 

300 CONTINUE 

301 RAO ( INPUT) ,LOA(NMAX) , INA (-1 ) ,AJP1 ( 22 0 ) , SL J ( 2 77 ) 

220 ENA(2) »STA(JA) , ENA ( 0 ) , ST A ( NH IT ) , LOA( 0.0 ) , ST A ( ELCST ) 
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BOX 19 



00 451 K= 1 , NMAX 
J=KH I T ( K ) 

DO 351 1*1,4 
351 SAVE( I,K)*B(I+3, J) 

00 43C 1=1,7 
430 COL D( I ) =0.0 
221 J=KH I T ( J A ) 

224 L0Q2 ( LB) ,LDL ( MASK 1 ) ,+AjP( L + 2 ) , EN A( 3) , -ST A (NPA IR ) 
00 245 K= 1 » KMAX 

LAC3( NUM) , INA2 ( 0 ) ,AJP(246) 

245 CONTINUE 

246 PAR=CSTANCE(K)*1 .0E-10 
DO 240 1=1,3 

B( I , J)=B( I , J) +B( 1+3, J)*T (K) 

LDA ( J A ) , I N A C - 2- ) , AJPI225) 

B( I » M ) *B ( I , M ) +B ( I + 3 ,M ) * ( T ( K ) -T I ME ) 

SL J ( 240) 

225 B( I ,M)=B( I,M)+B( I+3,M)*T(K) 

240 CONTINUE 

T I ME = T ( K ) 

LDA2ILB) ,SCL( MASK7) ,STA2 (LB) 

380 CALL INTERIM, J, PAR) 

TLAST( J)=TTIME-TSL+TIME 
805 DO 227 1=1,7 

COLD ( I ) =COLD ( I ) +DEV( 1,1) 

227 B ( I , J ) =DEV (1,2) 

COLD ( 7 ) =0.0 

EL0ST=B(7, J)-SAVE(4,JA)+ELCST 
SL J ( 22 1 ) 
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BOX 20 



T2=B(7,M)-ELOST 

AJP ( 356) ,AJP2(229) , LDA (0.0) , STA ( TEMP ) , S T A( T EM ) 

356 DO 357 1=4,7 

357 B( I,M)=0.0 

LDA ( T2 ) , A JP ( 229 ) 

DO 358 K= 1 , NMAX 
J=KH IT ( K ) 

520 DO 352 1=1,3 
TEMP=TEMP+SAVE( I , K ) »S AVE ( I ,K) 

TEM=TEM+B( 1+3, J) »e( 1+3, J) 

352 CONTINUE 

358 CONTINUE 
TEMP=TEMP/TEM 
STUPID=SQRTF( ABSF(TEMP) ) 

DO 355 K=2 , NMAX 

J=KH I T ( K ) 

521 DO 354 1=4,6 

354 B( I , J)=STUPIO*B( I,J) 

B ( 7 , J)=B(7,J)»TEMP 

355 CONTINUE 

229 TEM=NHI T 
TEM=1 .O/TEM 
DO 359 1=1,3 

359 B(I»M)=COLD(I)»TEM 

LDA ( T2 ) , A JP (231 ) , A JP3 ( 2 3 1 ) , LDA ( S AVE+ 4 ) , FSB ( ETHRE SH ) , A JP2 ( 503 ) 
LAC ( SAVE + 4) , FMU ( TURN ) , FAC ( ELOS.T ) , A JP ( L + 2 ) , A JP2 ( L+ 1 ) , SL J ( 50 3 ) 
LDQ51LB) , LOL ( MASK2 ) ,ARS( 27) , STA( IL ) , LDA( 0.0 ) , STA ( TEM) 

CALL SITE (IL) 

DO 501 1=1,3 
TEMP=I B ( I ) 

SCALE( I )=TEMP*A-B( I .M) 

FMU 1 ( SCALE ),FAO( TEM), STA (TEM) 

501 CONTINUE 

LDA(TEM) , FSB ( ALFA), AJP(5C3), AJP3( 503) 

TEM=SCRTF(ABSF(T2*C1 O/TEM) ) 

DO 502 1=1,3 

502 B( I+3,M)=TEM*SCALE( I) 

B ( 7 , M ) =T2 

GO TO 231 

503 LDA(NHIT) »INA(-1 ),AJP1 (505) 

DC 504 1=4,7- 

504 8 ( I ,M)=DEV( I, 1 ) 

GO TO 231 

505 TEMP=0 .0 

DO 329 1=4,6 

LDA 1 (COLD)fFMUI(COLD) , FAC ( TEMP ), STA ( TEMP) 

329 CONTINUE 
B ( 7 , M ) = T2 

TEMP=SQRTF(ABSF(C10»B(7, MI/TEMP) ) 

230 B( I ,M)=COLD( I )»TEMP 
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BCX 21 

231 ENA (0), ST A (NX), STM NY) , + LDA(NHIT) , AJ P 1 (L + 2) , STA( IN3) , SLJ1277) 
TLAST(M)=TTIME-TSL+TIME+TSL*PERCENT 
LOC ATE = 0 
TEMP=B( 7, M) 

FSB( ETHRESH) , A JP ( L+9 ) , + A JP3 < L + 8 ) , LOO ( M ASK2 ) , LCL5 ( LB ) , AJP (L + 4 ) 
+ARS(3),STA( I TEMP) ,LIU1 ( ITEMP) ,LDA1( LB) ,SCL (LB) ,STA1 (LB) ,LDA5(L8 ) 
SCL( MASK2 ) ,SST( 26) , STA5( LB) , +SLJ (514 ) .+LDA( TEMP ), FSB ( THERMAL ) 

A JP2 ( L+1 ) , SLJ ( L + 3) ,LDG(MASK2 ) ,LDL5(LB) , A JP 1 ( 5 1 3 ) , SL J ( L+2 3 ) 
+LDA5(LB),+SCL(2B),STA5( L 8 ) , LOC ( M ASK 2 ) , S TL ( I TEMP) • 

LDQ ( MASK 1 ) ,STL( J TEMP) ,LOA( JTEMP) , AJP 1 (L+6),+LDA( I TEMP) ,AJP( L+4 ) 
ARS ( 27 ) , STA( I TEMP ) , LI L 1 ( I TEMP ) ,L0A1( LB ) , SCL ( 4 B ) , ST A 1 ( LB ) , SL J ( L + 4 ) 
+ ARS( 15),STA( JTEMP) ,LIL1 ( JTEMP ),LDA1 ( LP. ) , SCL (MASK1 ),STA1(LB) 

LD A5 ( LB ) » SCL ( MASK3 ) , LDQ ( M A SK6 ) , STL ( I T EMP ) . SC L ( MASK 6 ) , ST A 5 ( L P ) 

LD A ( I TEMP) , A JP 1 ( L + 4 ) , LI LI ( NZ ) , ENA5(0 ) , STAl { LATER ) , IN I 1 ( 1 ) 

+ SIL1 (NZ) » +LDA ( T EMP ) , A JP ( L + 2 ) , AJ P2 (L + 1 ) , R AO ( M I ST AK E + 3 ) , E NA ( C ) 

STA ( LOCATE ) 

600 SIL51MJ) ,ENA ( 1 ), STA ( MOVE ), ENA ( 0 ), STA ( NOW ) 

CALL INVAC ( M J ) 

ENA(C) , STA (MOVE) , LO A5 ( LB ) , SC L ( MA SK7) , ST A5 ( L P ) ,L CC ( NO W ) , Q J P ( 5 1 C ) 
QLS ( 27 ) » STQ ( NOW 1 ) 

509 LDA5( LB ) ,SSU(MASK2) , STA5 ( LB ) , L I L 1 ( NOW ) ,LD01 ( LB) , LDL ( 48 ) , AJP 1 ( L + 3 ) 
-L0A1 (LB) »SST(4B)»STA1 (LB) ,SLJ(513), +ENA(0) , ST A ( MANY ) , ST A( NANY ) 

00 507 I = l , L M AX 

LOOK LB) ,LOL (MASK2 ) , SUB ( NOW 1 ) » AJP 1 (5 C7 ) ,LOL( MASKl ) , AJP( 508) 

ARS ( 1 5 ) » STA ( MANY ), LOO (NANY) ,QJP1 (507 ), STAl NANY) 

507 CONTINUE 

508 ENA1 (0) , AJP1 ( L+4) . + LIL1 ( NOW ) , LCA 1 ( LB ) , SCL ( 4 B ) , ST A 1 ( L6),LCC(NCWl ) 
SLJ15C9) ,LDA( MANY), SUP (NANY) , + AJP1 (L + 1 ) , RAO ( M I S T AK E+ 5 ) ,LDA ( NANY ) 
ALS ( 15) , STA ( I TEMP) ,ENA5(C) ,ALS( 15) ,STA( JTEMP ) , LO A5 ( LB ) , LOO ( M ASK 1 ) 
SSU( I TEMP), STA5(LB) , + LILl (NANY) , +LOA 1 ( LB ) , SSU ( JT EMP ) , ST A 1 ( L B ) 

ENA 1(0), SUB (MANY) ,AJP(51 3 ) , L I L 1 ( MANY ) , SL J ( L- 3 ) 

510 DO 51 1 K = 1 , JPMAX 

L0A31M00N) ,STA( I L) , LD A ( 0 . C ) , ST A3 ( OST ANCE ) 

CALL SITE ( I L ) 

DO 511 1=1,3 
TEMP= I B ( I ) 

TEMP=B( I ,M)-TEMP*A 

511 DSTANCE(K)=DSTANCE (K)+ TEMP* TEMP 

1 = 1 

DO 512 K=2, JPMAX 

L0A3 (DSTANCE ) ,FSB1 ( OST ANCE ) , AJP3 ( L+1 ) , SL J ( 5 1 2 ) , SIL3 ( IT ) , L I L 1 ( I T ) 

512 CONTINUE 

L0Q1 (MOON) ,STQ( NOW) ,0L S ( 27 ) , STC ( NCWl ) . SL J ( 509 ) 

LOA ( LOCATE ) , A JP 1 ( 5 1 6 ) , LD A ( T2 ) , A JP ( L+ 2 ) , A JP2 ( L+ 1 ) , R AO ( M I S T AK E + 3 ) 
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BOX 22 



514 ENI 1 (2) 

515 LIJA1 jKHIT ) ,STA( IT) , LIL2( IT) 

FSB? ETHRESH) , AJP (L + 9) , +A JP 3 ( L*8) ,LDQ (MASK2 ) .LCA2 (LB) , STL ( I TEMP ) 

SCL ( MASK2 ) »SST(2B) , STA2 ( LB ) , LOA ( ITEMP),AJP(233),ARS(3),LIU3IITEMP) 
LDA3 ( LB ) * SCL ( 4B ) , ST A3 ( LB ) , -*-S L J ( 2 3 3 ) , +LCA( T1 ) , FSB (THERMAL ) 

AJP(L+1) ,AJP2(L+17) ,LDC(MASK2) ,LDA2( LB) , SCL ( 2fi ) , ST A2 ( LB ) 

STL ( I TEMP) ,LDG(MA$K1 ) t STL ( J TEMP) , LDA ( J T EMP ) , A JP ( L +8 ) , ARS (15) 

ST A ( JTEMP ) , L I L3 ( JTEMP ) ,LCA3( LB) , SCL( MASK 1 ) , STA3 ( LB ) , LDA2 ( LB ) 
SCL(MASK3),SST(MASK6) , STA21LB) , L I L3( N Z ) , ENA 2 ( 0 ) , ST A3 ( L AT ER ) 

RAO(NZ) • SL J ( 2 33 ) , LDA ( I TEMP ) , ARS ( 27 ) , STA ( I TEMP ) ,L I L 3 ( I TEMP ) 

LDA3(LB) ,SCL(4B) , + SLJ (L-E) ,♦ LOG (MASK 2) , LDA2 ( LB ) , SST ( 28 ) , ST A2 ( LB ) 
LOL2 (LB) ,AJP1 (233) ,SIL5( LOCATE) ,SIL2 ( N A V Y ) , L I L5 ( NAVY ) ,SIL1 (MAO) 
SLJ( 600) 

516 LIL2(NAVY) ,LI L5( LOCATE) , LIL1 (MAD) 

233 ENA 1 (0) ,SUB(NMAX) ,A JP ( 234 ) , I NI 1 ( 1 ) ,S L J ( 5 15 ) 

BOX 23 

234 LOA ( NHIT) , AJP1 (241 ) , LDA( C.O) , STA ( ZZ+ 6 ) , SL J ( 283) 

241 00 282 NN= 1 « NMA X 

J=KH I T ( NN ) 

DO 280 K=1,KMAX 

LAC 3 ( NUM ) » INA2( C ) »A JP ( 28 1 ) 

280 CONTINUE 

281 LDA(NN) , I NA ( — 1 ) , + A JP ( L+2 ) , LO A3 ( T ) , -S T A ( T I ME ) 

DO 235 1=1,3 

235 B( I , J)=B( I , J)-( TIME + ZZ( 6 )*1 ,0E+14)*B ( 1+3, J) 

282 CONTINUE 

BOX 24 

283 GO TO (277,277,236,277) ,NPAIR 

236 DO 238 K=2,NMAX 
J=KH I T ( K ) 

T2=B(7, J) 

FSB (ETHRESH) , A JP ( 52 6) , A J P3 ( 5 26 ) , LDG2 ( LB ) , LDL ( M A SK2 ) , AJ P( 2 3 8 ) 

STA ( IL) , LDL (MA SKI ) , AJP ( 2 38 ) , STA ( LL ) , LDA2 ( LB ) , SC L ( MASK3 ) , STA 2 ( LB ) 
LDA(LL) ,ARS( 1 5) ,STA(LL) 

T4=B(7,LL ) 

FSB (ETHRESH) ,AJP3(L + 7) , L I L 1 ( LL ) , LDA1 ( LE ) , SCL ( MA SK3 ) , S T A 1 ( LB ) 

LDA ( IL), ARS( 3) , S T A ( I L ) , L IU 1 ( I L ) , LD A1 ( L B ) , SCL ( 48 ) , S T A 1 ( LB ) , SL J ( 23 e ) 
LIL 1 (LL) ,LDA1 (LB) ,SCL (M.ASK1 ) ,STA1 (LB ) , LDQ( MASK7 ) , STL ( IT) ,LDA( IT) 
AJP (238) ,LDA( IL) , ARS (27) , + STA( IL) 

CALL SITE ( IL) 

DO 237 1=1,3 
TEMP= I B ( I ) 

237 B ( I , LL ) =A*TEMP 
GO TO 238 

526 LDQ2 (LB) ,LDL( MASK1 ) ,AJP( 238) ,ARS( 15) , ST A ( LL ) , L I L 1 ( LL ) 

T4=B( 7, LL ) 

FSB ( ETHRESH) , AJP 3 (238) , LCA1 ( LB ) , SCL( MASK3 ) , STA 1 (LB) ,LDA2 (LB ) 

SC L ( MASK 1 ) ,STA2 (LE) 

238 CONTINUE 

BOX 25 

277 LDA ( IN3) ,AJP1 (200) , LDA ( I NMIN) , SUB( IN MAX) , AJP (302 ) , AJP2( 30 3 ) 

302 M=LAST( INMIN) 

RAO ( INMIN)„LDQ5(LB) ,LDL(MASK6) ,AJPl( 277) ,SLJ(20C) 
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BOX 26 



303 



304 

525 

350 

305 



261 

262 
2 50 

251 

252 

254 

255 
307 



331 

7998 

8000 



t 309 
310 



31 1 
312 



314 

315 
257 



332 

320 

316 



9000 



DO 3C5 J = 1 » I NM AX 
K = L AST ( J ) 

LDQ3 ( LB ) , LDL ( MASK6 ) , AJP1 (52b ) 

DO 304 1=' 3 
B ( I f K ) =B ( 

LDAi ( LB ) f 
I F ( ENERGY 
ENERGY=B ( 

I EN ERGY = K 
CONTINUE 

BOX 27 

LDA(NZ), INA(-l) » A JP ( 262 ) ,AJP3(262),STA(IT1) 

DO 261 I T= 1 « I T 1 
M J =L A T ER ( IT) 

CALL INVAC(MJ) 

ENA(O) ,STA( IU)*SIU6(IU) LDA(MJ) RAD(IU) 

WRITE TAPE 6»IU»(B(I*MJ) 1 = 1 , 3 ) 

CONTINUE 

' BOX 28 

DO 7998 J = 1 , INMAX 

LDA2 (LAST) ,STA( I TEMP) , LIL5( l TEMP ) , LD Q5 ( L B ) ,LDL( MASK7 ) , AJP1 ( 7999) 
DO 252 1=1,3 

IF ( B ( I ,M) -AR ( I ) ) 251.251,254 
I F ( B ( I , M ) ) 255,252,252 
CONTINUE 
GO TO 307 

ENA 1(0),ALS(1 )»STA(NFACE)*SLJ(252) 

ENAl ( 0 ) , ALS ( 1 ) , I NA (- 1 ),STA(NFACE) ,SLJ(252) 

LDA( NFACE) ,AJP( 7998 ) « ALS ( 36 ) , INA6(0) ,ALS(3),STA(ITEMP),LDA5(LB) 
LDQ ( MASK4 ) ,SSU( I TEMP ) , LDC ( MA SK2 ) , STL ( I TEMP) , LDQ( MASK 1 ) ,STL ( JTEMP ) 
SCL (MASK3 ) ,STA5 (LB) ,ENA( C) ,STA( NFACE ) , LDA( JTEMP) , AJP 1 ( 33 1 ) 

LDA( I TEMP) ,ARS(27),STA( I TEMP) ,LIL1 (I TEMP) ,LOAl (LB) ,SCL(4B) 

STA1 (LB) , SL J ( 7998 ) 

ARS( 15) ,STA( JTEMP) , L I L 1 ( JTEMP) ,LDAl ( LB) , SCL ( MASK 1 ) , STAl ( LB) 

CONTINUE 

CONTINUE 

THIS IS THE END OF THE MAJOR DO LOOP ON THE PARTICLES 
BOX 29 

ENA ( 0 ) , ST A ( IT ) ,SIU6( IT) 

DC 9C00 L = 1 , LM AX 

LDQ4 (LB) ,LDL( 2B) , AJP (309 ) ,SIL4( IT) 

WRITE TAPE 6,IT,(R(I,L) 1=1,3) 

LDQ4(LB) ,LDA(MASK8) ,STL4(LB) ,LDL(77771B) ,AJP1 (9000) 

LDL ( MASK7) ,AJPl ( 320 ) 

DO 310 1=1,3 

B(I,L)=B(I,L)+B(I+3,L)*TSL 
DO 312 1=1,3 

IF ( B( I , L ) — AB ( I ) ) 311,311,314 
I F ( B ( I , L ) ) 315,312,312 
CONTINUE 
GO TO 257 

ENAl ( 0) ,ALS( 1 ) ,STA( NFACE ) ,SLJ(312) 

ENAl (0) ,ALS( 1 ) , I NA ( — 1 ) , ST A (NFACE ) , SL J ( 3 1 2 ) 

LDA ( NFACE ), AJP (320) , ALS ( 36 ) , I NA6 ( 0 ) , ALS ( 3 ) , S TA ( I TEMP ) , LD A4 ( LR ) 

LDQ ( MASK4 ) , SSU ( I TEMP) , LDC ( MASK 2 ) , STL ( I TEMP) ,LDQ(MASK1),STL( JTEMP) 
SCL ( MASK 3 ) ,STA4(LB) ,cNA(C) , STA (NFACE ) ,LDA( JTEMP ) ,AJPK 332) 

LDA( I TEMP ) , ARS127) ,STA( I TEMP) ,LIL1 (I TEMP) ,LCA1 (LB) ,SCL(4B ) 

STAl (LB) , SL J ( 320 )' 

ARS ( 15) , STA ( JTEMP) ,LIL1 ( JTEMP.) ,LCA1( LB), SCL ( MASK 1 ) .STAl ( LB) 

I F ( ENERGY-B ( 7 , L ) ) 3 1 6 , 90C C , 9000 
ENE RG Y=B ( 7, L ) 

I ENERGY=L 
CONTINUE 

THIS IS THE ENC OF THE MAJOR DC LOOP ON THE TIME STEPS 



-B ( 7 , K ) ) 350,305, 



7 , K ) 



305 
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BOX 30 



F 1 3. 7,7H JIFFYS ) 
F13.7,7H JIFFYS) 



8001 WRITE TAPE 8 , I X , I Y , I Z , AB , N'L I NE , N PL AN E , MAXL A T , LM A X , A , S PH I , FM A S S , 

1 E THRESH, THER MAL, EPB, LULL, ALPHA, 8 ETA, N ,MNTS»TTIME»TSLI»NTSC» 
20ENERGY,TENERGY, MISTAKE, YENTRY, ZENTRY ,N A , EMAX, EFOUND 

WRITE TAPE 8, ( LB ( L ) , (R( I ,L) 1=1,7) L=1,LMAX) 

REWIND 8 

♦ SL S 1 ( L + 1 ) ,+SLJl (L + l ) , SL J ( L + 2 ) ,SLJ2( L + 1 ) ,SLJ3(999B) 

WRITE OUTPUT TAPE 3, 9951,TIIME 

9951 FOR MAT ( 1 H 1 , 1 9H THE TOTAL TIME IS 
WRITE OUTPUT TAPE 3, 9952, TSL 

9952 FORMAT ( 25H THE TIME STEP LENGTH IS 

I F ( N ) 9974,9974,9973 

9973 WRITE OUTPUT TAPE 3 ,99 53 , ENERGY , N 
GO TO 9972 

9974 WRITE OUTPUT TAPE 3 , 9953 , ENE RGY , MNTS 

99530F0RMATI9H ENERGY= F15.7,15H ELECTRON VOLTS/14, 27H TIME STEPS HAV 
IE BEEN TAKEN ) 

9972 WRITE OUTPUT TAPE 3 , 9954 , T EN ERGY , OEN ERGY 

99540 FORMAT! 2 1H THE TOTAL ENERGY IS F12.7.15H ELECTRON VCLTS/24H THE C 
1RIGINAL ENERGY IS F12.7.15H ELECTRON VCLTS) 

I F ( M I STAKE ( 1 ) ) 9955,9955,9956 

9955 WRITE OUTPUT TAPE 3,9957 

9957 FORMAT ( 30H THE ENERGY CHECK DID NOT FAIL ) 

GO TO 9965 

9956 WRITE OUTPUT TAPE 3 , 99 59 , M I ST AK E ( 1 ) 

9959 FORMAT ( 24H THE ENERGY CHECK FAILED I4,16H DIFFERENT TIMES) 

9965 IF ( MI STAKE ( 2 ) ) 9961,9961 ,9962 

9961 WRITE OUTPUT TAPE 3,9963 

9963 FORM AT ( 62H NO INTERSTITIALS OR VACANCIES 
1C0RNER ) 

GO TO 9805 

9962 WRITE OUTPUT TAPE 3 , 9960 , M I S T AKE ( 2 ) 

9960 FORM AT ( 60H AN INTERSTITIAL CR VACANCY CCCURRED 
1RNER 14, 1 6H DIFFERENT TIMES) 

9805 I F { M I STAKE ( 3 ) ) 9800,9800,980 1 

9800 WRITE OUTPUT TAPE 3,9803 

9803 FORMAT ( 36H NO NEGATIVE ENERGIES WERE COMPUTED 
GO TO 107 

9801 WRITE OUTPUT TAPE- 3 , 9804 , M I S T AKE ( 3 ) 

9804 FORMAT ( 3 1H A NEGATIVE ENERGY WAS COMPUTEC 14, 

1 /) 

107 I F ( M I STAKE ( 4 ) ) 170,170, 172 

170 WRITE OUTPUT TAPE' 3,171 

171 FORMAT ( 37H NO FALSE INTERSTITIALS WERE FORMED 
GO TO 174 . 

172 WRITE OUTPUT TAPE 3 , 1 7 3 , M I ST AKE ( 4 ) 

173 FORMAT (13, 34H FALSE INTERSTITIALS WERE FORMED /) 

174 WRITE OUTPUT TAPE 3 , 1 75 , MI ST AK E ( 5 ) 

175 FORMAT ( 10H WE FORMED ,IL,22H TRIPLE INTERSTITIALS ) 

WRITE OUTPUT TAPE 3 , 1 08 , I X , 1 Y , I Z 

v 1080F0RMAT(35H THE DIMENSIONS OF THE LATTICE ARE 



OCCURED ALONG AN ECGE OR 



ALONG AN ECGE CR CO 



/ ) 



16H DIFFERENT TIMES 



/) 



106 

X 



1 12, 7H UNITS 
212, 7H UNITS 
WRITE OUTPUT 
106 FORM AT ( 87H 
1 Ml 

DO 1025 L = 

1025 WRITE OUTPUT 
103 FORMAT! IX, I 4, 7F10. 4, 01 9) 

9998 ENA610),AJP(9990),SLS1( 1 

9990 WRITE OUTPUT TAPE 3,9991 

9991 FORMAT!/// 32H THE RUN WENT 

9999 JT=- 1 

WRITE TAPE 6 , JT 

REWIND 6 

END 



/ 1 3H Y DIRECTION 12, 7H UNITS 
) 

TAPE 3, 

L 

E LB ( L ) / ) 

1 , LMAX 

TAPE 3 , 1 0 3 » L 



// 1 3H X DIRECTION 
/ 1 3H Z DIRECTION 



VX 



V Y 



( 6 ( I , L ) 1= 1 ,7) ,OCTAL(L ) 



16) 



THE FULL 0 1 STANCE 
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GOX 31 



SUBROUTINE SITE(V) 

DIMENSION R ( 7 « 1 5CC ) f LB ( 1 5CC ) .OCTAL (1 5CC),FI(3),IB!3),T(20) 
DIMENSION TMINI2C) , I FACE (3 ), MI STAKE! 1C ) ,MOCN( 15 ) ,NGON( 101.MH1 
DIMENSION XC ( 6 ) , XP ( 20 ) , Z Z ( 1 5 ) , CON ! 4 ) , C E V ( 7 , 2 ) , A B ( 3 ) 

DIMENSION DST ANC E ( 20 ) 

COMMON B.LB.OCTAL.FI, I B , T , TM I N , I FACE ,M I S T AK E , MOON .NOON , W , XC , X P , Z Z 
COMMON CON.DEV, AB, A.NPLANC ,N LINE, LMAX,MASK1,MASK2,MASK3, MAS K4 
COMMGN I TEMP, JTE MP, MT EMP ,N TE MP , IN2,IN3,IL,JA,JB,JC,JK,JL,JM,JP,JR 
COMMON JT,JW,MM,KMIN,MAX,TEMP,NI,JPMAX,NIF,IX.IY,IZ,NN,MAXLAT,JJ 
COMMON IP,TE,TEM,C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,C12,T1,T2 
COMMON DSTANCE, FMASS, 13, TSL 
COMMON MASK5,MASK6,MASK7, MAS K8, MOVE, NCV* 

EQUIVALENCE ( XP , TM I N ) , { L e , CC T AL ) 

JT = M 

LDA( JT), I N A ( — 1 ) , ENG ( OR ) ,DVI (NPLANE ) , STA(IR + 3),LLS!48),ENC(0) 

DVI (NLINE) ,STA! IB + 2 ) ,CLS ( 1 ) ,STQ( IB + 1 ) , ADD ( I e + 3 ) , ENQ! 1 ) , STL ( JT ) 
LDA( JT) , AJP( 234C ) , R AO ( IB+1 ) 

2340 RETURN 
END 

BOX 32 

SUBROUTINE INVAC(MJ) 

DIMENSION B(7, 15C0) ,LB(1500) .OCTAL (1 500) ,FI (3) , IB! 3) ,T(20> 
DIMENSION TMIN! 20) , I FACE (3) .MISTAKE! 1C) .MOON! 15) .NOON! 10 ) ,W(4 ) 
DIMENSION XC ( 6 ) , XP ( 2 C ) , Z Z ( 1 5 ) ,CCN ( 4 ) , DE V ! 7 , 2 ) , AB ( 3 ) 

DIMENSION DSTANCE ( 20 ) 

COMMON B, LB, OCTAL, FI , I R , T , TM I N , I FACE , M I S T AK E , MOON , NOON, W , XC.XP.ZZ 
COMMON CON, DEV, AB, A, NPLANE .NLINE, LMAX.M A SKI , M ASK2 , M A SK 3 , M A S K 4 
COMMON ITEMP, J TEMP, MT EMP ,N TEMP, I N2, I N3, IL,JA,JB,JC,JK, JL.JM.JP, JR 
COMMON JT, JW,MM,KMIN,MAX,TEMP,M ,JPMAX,NIF,IX,IY,IZ,NN,MAXLAT,JJ 
COMMON IP,TE,TEM,C1 ,C2,C2,C4,C5,C6,C7,C8,C9,C10,C1 1 ,C12,T1,T2 
COMMON DSTANCE, FMASS, 13, TSL 
COMMON MASK5.MASK6, MASK7 .MASK8.M0VE, NOW 
EQUIVALENCE ! XP , TMI N) , ! LB, OCT AL ) 

2360 IN2=0 
MM= M J 

LILKMM) ,LDA1 (LB),SCL (MASK2) ,SCL IMASK6) , SST ( MASK7 ) ', ST A 1 ( LB ) 

2361 NIF = 0 

LDA {MIST AKE + 2 ) ,STA(MISTAKE+10) 

DO 2370 1=1,3 
I F AC E ! I ) =0 

I F ( AB { I ) — A-B ( I. KM) ) 2375,2375,2365 
r 2375 ENA1 !CB) , ALS! 1 ) , ST A 1 { I F ACE ) , RAO ( N IF ) ,SLJ!2370) 

2365 IF! B! I ,MM) — A) 2380,2370,2370 

2380 ENA1 (OB) ,ALS< 1 ) , INA ( - 1 ) ,STA1 {I FACE), RAC (N IF ) 

2370 CONTINUE 

DO 2387 1=1,3 
IB! I )=8( I ,MM) /A 
TEMP=I B ! I ) 

I F ( B ( I , MM ) / A— TEM P—0 . 5 ) 2 367,2386,238 6 

2386 RA01 ( IB ) 

2387 CONTINUE 
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BOX 33 



IL=IB(3)*NPLANE+IB(2)*NLINE+IB(1)/2+l 

LDA( IB+1 ) » ADD ( IP, +2) ,ADD( I P + 3) ,ENG( 1) ,STL(ITEMP) ,LDA( ITE”P) 
AJP1(2410),STQ(JM),SLJ(2415) 

2410 JK=0 

2415 LDA ( I B + 2 ) » ADD ( IB+3),ENC(1),STL(JL),RSC(JL),LCA(JM) ,ADD( JL) 

+ AJP3 { L + 3 ) * LDA( IL ) * STA(NCW) , LDA (MOVE ) , A JP 1 ( 274 1 ) , SLJ ( 24 3 2 ) , RAC ( I L ) 
2432 J A = C 

DD 2115 1*1 » 3 

LDAl ( I FACE) ,AJP( 21 15) , AJP3I2 11 51 ,RA0 ( JA) , INA(-1 ) ,STA( JW) 

AJP1 ( 2120) ,LDA1 ( I FACE) , S TA ( JB ) , ENA ( 1 ) , S TA ( JD ) 

2115 CONTINUE 

LDA( JA) , AJP( 2125) ,SLJ(2140) 

2120 RAOINISTAKE+2) »SLJ(214C) 

2125 ENA ( 0 ) , STA ( JB ),STA(JD),LCA(JM),STA(JC),SLJ(214l) 

2140 JC= 1 

214 1 ENA ( 1 ) ,STA( JK ) ,LDA ( JM) ,MLI ( 6B) ,MUI ( JD ) , I NA( 1 ) , ADD( JB) , ADC ( JC ) 

STA(JW),INA(-2),AJP(2447),AJP2(2448), ENA (6) , STA ( JPMA X ) , SL J ( 2 700 ) 
2447 ENA{ 1 5B) , STA (JPMAX) , SLJ ( 2700 ) 

24 48 LDA ( JW) , INA(-1 IB ),AJP3 (2449) , ENA ( 1 IB ), STA( JPMAX) ,SLJ (2700 
2449 ENA(5),STA( JPMAX) 

2700 CO 2698 1=1,10 

2698 NOON ( I )=0 

DO 2699 1=1,15 

2699 M 0 0 N ( I )=0 

BOX 34 

ENI 1 ( 1 ) 

DO 2475 JP=1 , JPMAX 

GO TO (2701,2702,2703,2704,2705,2706,2707,2708,2709,2710,2711, 
12712,2713,2714). JW 

2701 GO TO (2455,2451,2452,2453,2454,2465) ,JP 

2702 GO TO (2455,2451,2470*2456,2452,2453,2457,2458,2459,2454,2460, 
12461 ,2462), JP 

C2500 THIS SECTION IS USED WHEN THE PARTICLE SLCWS DOWN WITHIN CNE UNIT 
OF AN EDGE. WE CAN NOT CONSIDER ALL THE NEAREST NEIGHBORS IN A 
CASE LIKE THIS, ONLY CERTAIN SELECTED ONES, DEPENDING ON WHICH 
FACE IT IS THAT WE ARE NEAR. 

2703 GO TO (2455,2451 ,2452,2453,2454) ,JP 

2704 GO TO (2451 ,2452,2453,2454,2465) ,JP 

2705 GO TO (2455,245 1 ,2453,2454,2465) ,JP 

2706 GO TO (2455,2452 ,2453,2454 ,2465) ,JP 

2707 GO TO (2455,2451 ,2452,2453,2465) ,JP 

2708 GO TO (2455, 2451, 2452, 2454, 2465), JP 

2709 GO TO (2455,2470,2456,2457,2458,2459 ,2460,2461 ,2462) ,JP 

2710 GO TO (2455,2451 ,2452, 2453, 2454, 2458, 2459, 2461, 2462), JP 

2711 GO TO (2455,2451 ,2470, 2453, 2457, 2453, 2454, 2460, 2461), JP 

2712 GO TO (2455,2452,2453,2454,2456,2457 ,2459,2460,2462) ,JP 

2713 GO TO (2455,2451 ,2452,2453,2456,2457 ,2458,2459,2470) ,JP 

2714 GO TO (2455,2451 ,2470,2456,2452,2454 ,2460, 2461, 2462), JP 
2451 LDA( JL) ,ADD INLINE ) , ADD ( I L ) ,SLJ(2435) 

24 52 LDAl IL) , SUB INLINE) ,ADD( JL) ,SLJ(2435) 

245 3 LDAl IL) , ADD( NPLANE) , ADD ( JL ) , SL J ( 24 35 ) 

2454 LDA ( I L ), SUB ( NPLANE ), ADC ( JL ), SLJ ( 2435 ) 

2455 LDAl I L ) , SLJ ( 2435 ) 

24 56 LDA ( IL) ,SUB( NLINE) . INA( 1 B) , ADD( JL) ,SLJ(2435 ) 

2457 LDAl IL) , ADD (NPLANE ) ,INA( IB) , ADD( JL) , SLJ (2435 ) 

2458 LDAl IL), ADD (NPLANE) ,ADD(NLINE) ,SLJ(2435) 

24 59 LDAl I L ) ,ADD ( NPLANE ) , SUB ( NL I NE ) , S L J ( 2 435 ) 

2460 LDAl IL) , SUB (NPLANE) ,INA( 1 B ), ADC ( JL ), SLJ ( 2435 ) 

2461 LDA ( IL) ,SUB( NPLANE) , ADC( NL I NE ) , SLJ (2 4 35 ) 

2462 LDAl IL) , SUB ( NPL ANE ) , SUB ( NL IN E ) , S L J ( 2 435 ) 

2465 LDAl ID, INA(-l) ,SLJ(2435) 

2470 LDAl IL) , ADD ( NL I N E ) , ADO ( J L ) , I NA ( 1 ) 

2435 STA( I TEMP) ,AJP (2475 ), AJP3 (24 75) ,SUB( MAXLAT) , +AJPI L+l ) ,AJP2(2475) 
LDAl I TEMP ) ,STA1 ( MOON) , INI 1 ( 1 ) 

2475 CONTINUE 
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BOX 35 



INI 1 ( - 1 ) , SILK JPMAX) , LD A ( MC V E ) ,AJP1( 274 1 ) 

2480 EN A ( 1 ) , S T A ( J P ) 

00 2485 K=l, JPMAX 

LDA3 ( MOON) , STA ( I TEMP ) , L I L 1 ( I TEMP ) , LO C 1 ( L B ) , LDL ( 4 B ) , A J P 1 ( 24 8 5 ) 

ENQ 1 (0) * LIL1 ( JP) , STQ 1 ( NOCK) ,RAQ( JP) 

2485 CONTINUE 

LDA( JP) , INA(-1)»STA(JK) 

BOX 36 

LDA ( JK ) , INA(-1 ) , AJP( 2490) , AJP2 (2510) , LC A ( I N2 ) , A J P 1 ( 2720 ), SLJ ( 2600 ) 
2490 LDA (NOON+1 ) , STA( IL) 

2492 CALL SITE ( ID 

2493 DO 2495 1=1,3 
T EM P = 1 8 ( I ) 

B ( I ,MM ) =A*TEMP 
2495 CONTINUE 

ENA(3) ,STA(NI ) ,L I LI ( IL) , LIL3(MM) ,LOA ( IL) ,ALS(27) ,STA( I TEMP) 

LDA3 ( LB ) ,LDQ( MA SK2 ) , S SU ( IT EM P ) , ST A3( LB ), L0A1 (LB),SST(4B),STA1(LB) 
LDA { IN2), AJP1 (2800 ) 

RETURN 

BOX 37 

2510 DO 2512 K=1 , JK 

LDA3 ( NOON ) ,ST A ( I TEMP ) , LD A ( C . 0 ) , STA3{ CST ANCE ) 

CALL SITE ( I T EMP) 

DO 2512 1=1,3 
TE MP= I B ( I ) 

TEMP=TEMP*A 

2512 DSTANCE ( K ) =DSTANCE ( K ) + ( TEN P-B ( I , MM) ) * { TEMP-B ( I , MM ) ) 

EN I 2 ( 1 ) 

DO 2514 K = 1 , J K 

LDA3 (DSTANCE ) , F S B2 ( DST ANCE ) , A J P ( 2 5 1 4 ) , A J P2 ( 2514),SIL3(ITEMP) 

LIL2( ITEMP) 

2514 CONTINUE 

EN I 5 ( 1 ) ,LDA2(N00N) ,STA5( MOCN ) , INI 5(1 ) 

DO 2516 K= 1 , JK 

LDA2( DSTANCE) , F SB3 ( DST ANCE ) , + A JP ( L + l ) , A J P2 ( 2 5 1 6 ) , +S I L2 ( I T EMP ) 

' LAC( ITEMP) , I N A3 ( 0 ) , AJ P ( 2 5 1 6 ) ,LDA3( NOCN) , STAS ( MOON ) , I N I 5 ( 1 ) 

2516 CONTINUE 

I N I 5 ( - 1 ),ENA5(0) , IN A ( - 1 ) ,AJP1 (L+3) ,+ LDA2 ( NOON) , S T A ( I L ) , +SL J ( 2492 ) 
+ENA5( 0) , STA (MMAX) 

DO 2518 K=1 , MMAX 

LDA ( MM ) , SUB3 ( MOON ) ,AJP(25?0) 

2518 CONTINUE 

2520 ENA 3 ( 0 ) , AJ P { L + 3 ) ,LDA3 ( MOCN ), STA ( IL),+SLJ(24 92) , +ENQ6 ( 0 ) , LCL ( M ASK 5 ) 
+AJP1 (L+2) ,LDA(M0CN+2) , S TA ( I L ) , SL J ( 2 492 ) ,+LDA(MCCN+l ) , ST A ( I L ) 

SLJ ( 2492 ) 
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BCX 38 



2600 

2605 

2620 



1700 

2621 



2625 



2630 

2635 

2640 

2645 



2650 

2655 

2660 



2665 

2675 

2680 

2686 



T 1 = 0 . 0 

DO 2605 1=4,6 
T 1 =T 1 +8 ( I,MM)«B( I ,MM) 

ENA( 1 ) ,STA( JA) ,LDA(T1 ) ,AJP1 ( 2620) ,ENA(0),STA(JA),SLJ(262C) 
T2=0 . 0 



DO 2645 K= 1 » J PMA X 

LOA 3 ( MOON ) t STA( I TEMP) , ALS( 27 ) , STA( JT EMP ) 

DC 1700 1=1 , LM AX 

LD01 (LB) ,LDL(MASK2) ,SU8( JTEMP) ,AJP1 ( 17C0) ,LDL( MASK) ) , AJP (2621 ) 
LDA ( TSL 1 , FMU ( 1 . OE + 3 ) , STA3 ( TM I N ) , SL J( 264 5 ) 

CONTINUE 

CALL SITE(ITEMP) 

DO 2625 1=1 ,3 
T EMP= 1 8 ( I ) 

FI ( I ) = A*T EMP 
LDA( JA) , AJP 1 (2635) 

DO 2630 1=1,3 

TMI N(K) = TMIN( K) + (B( I ,MM)-F I ( I ) )* ( B ( I , MM) -FI ( I ) ) 

GO TC 2645 
DO 2640 1=1,3 

T2 = T2+(8( I.MMJ-F I( I ) )*B( 1 + 3, MM) 

TMIN(K)— T2/T1 
CONTINUE 

ENA ( 1 ) ,STA( KM IN) , LDA ( JPMAX ) , I NA ( -1 ) , ST A ( JK ) 

DO 2660 K= 1 , JK 

IF(TMIN(K)-TMIN(K+1 ) ) 266C ,2660,2650 
I F ( TMIN(K+1 J-TMIN(KMIN) ) 2655,2660,2660 



KM I N = K+ 1 
CONTINUE 

LI LI (KMIN) , LDA 1 (MOON) , ST A ( I P ) , LD A ( TS L ) , F MU ( 1 00 . 0 ) , S T A ( TE M ) 

LDA 1 (TMIN) ,+FSB( TEM) ,AJP3(L+2) ,ENI1 ( 0) , SL J { 2685 ) 

CALL SITE (IP) 

DO 2665 1=1 ,3 
TEMP= I 8 ( I ) 

FI ( I ) = A* TEMP 

T ( I ) = ABSF ( 8 ( I , M M ) - F I ( I ) 1 

LCA( JB) , INA( 1 ) . ARS( 1 ) ,STA( ITEMP) ,LIL 1 ( ITEMP) ,LDA(C.O) ,STA1 ( T ) 
ENA ( 1 ) .STA(MAX) 

DO 2675 1=1,2 

LOA 1 ( T) ,FSB1 ( T+ 1 ) , AJP (2 675 ) , AJP2(2-67 5) ,S ILl ( MAX ) ,RAO( MAX ) . 
CONTINUE 

DO 2680 I = 1 , LMAX 

LDQ1 (LB) , LDL ( MASK2 ) , ARS ( 27 ) , SUB ( I P ) , AJP (2685) 

CONTINUE 

STA(MM) , ENA ( 1 ) , S T A ( I N2 ) , S L J ( 236 1 ) 
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BOX 39 



2800 LDA( IP) *STA ( IL) « LDA( MJ ) , S T A ( MM ) , E N A( C ) , S TA { I N2 ) , S L J ( 249 2 ) 

THIS IS THE SECTION OF INVAC THAT MAKES THE INTERSTITIALS. THE REST OF 
THE SUBROUTINE IS CONCERNEC ONLY WITH VACANCIES 
2720 TEMP=B(MAX,MM)-B (MAX,HJ) 

CALL SI TE ( IP ) 

DO 2725 1=1,3 
T2= I B ( I ) 

B { I , MM ) = A*T2 
2725 B( I , M J ) = B ( I, MM) 

I F ( TEMP ) 2735,2730,2730 
2730 B(MAX,MM)=B(MAX,MM) +0.5*A 
B(MAX,MJ )=B (MAX, MJ)— 0.5*A 
GO TO 2740 

2735 B(MAX,MM)=B(MAX,MM)-0.5*A 
B(MAX,MJ )=B(MAX,MJ)+0.5*A 

2740 CONTINUE 

LDA { IP) ,ALS (27) , STA { I TEMP ) , LDA { MM ) , A L S < 1 5 ) , A CD ( I TEMP ) , STA( JA) 

LD0(MASK3) , L I L 3 ( M J ) , LD A3 ( LB ) , SSU ( J A) , ST A3 ( LB ) , ENA3 ( 0 ) , ALS ( 1 5 ) 

ADD ( I TEMP) ,STA( JA) , LI LI ( MM) , LDA 1 (LB),SSU(JA),STA1(LB),LIL1(IP) 

LDA1 (LB) ,SST(4E),STA1 (LB) ,LUA(MM),SUB(MJ ),AJP1 (2741 ) 

RAO ( MI STAKE+4 ) 

2741 RETURN 
END 



BOX 40 

SUBROUTINE I NTER ( MO , I I , S ) 

DIMENSION B(7, 1 500) ,LB{ 1500) ,0CTAL(1 500) ,FI (3) , I B< 3) ,T(20) 
DIMENSION TM IN ( 20 ) , IFACE (3) , M I STAKE! 10) ,MOON (15), NOON( 10 ) , K(4 ) 
DIMENSION XC ( 6 ) , XP ( 20.) , ZZ ( 1 5 ) , CCN ( 4 ) , DE V ( 7 , 2 ) , AB ( 3 ) 

DIMENSION DST ANC E ( 2 0 ) 

COMMON B, LB, OCTAL, FI ,IR, T , TM I N, I FACE ,M I STAKE , MOON, NOON ,W , XC , XP , Z Z 
COMMON CON, DEV, AB,A .NPLANE ,NLINE ,LMAX, MASKl , MASK2, MASK3 , MASK4 
COMMON I TEMP , JTE MP , MT E MP ,NTE MP , I N2 , 1 N3 , I L , J A , JB , JC , JK , JL , JM , JP , J R 
COMMON JT , JW, MM, KM I N, MAX , TEMP , NI ,JPMAX,NIF,IX.IY,IZ,NN.MAXLAT,JJ 
COMMON IP,TE,TEM,Cl,C2,C2,C4,C5,C6,C7,C8,C9,Cl0,Cn,Cl2,Tl,T2 
COMMON DSTANCE,FMASS, 13, TSL 
COMMON MASK5,MASK6,MASK7,MASK0,M0VE, NOW 
EQUIVALENCE ( XP , TM I N ) , ( L E , CC TA L ) 

MM = MC 

DO 1C 1=1,6 

XP ( I )=B( I , MM ) -B ( I, I I ) 

FMU ( 0 . 5 ) STA ( ZZ+2 ) . 

10 XC ( I ) =ZZ ( 2 ) 

ZZ ( 9)={ XP{4) *XP (4 )+XP( 5 )*XP(5)+XP(6) «XP(6))*1.0E8 
FMU ( C4 ) STA(ZZ+10). 

ZZ ( 1 )=SQRTF( ABSF(ZZ(9) ) ) 

ZZ ( 10) = 1 .O/ZZl 10) 

ZZ ( 1 1 )=S»S 
ZZ( 7 ) =C1 
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BOX 41 



15 ZZ(2)=C2/EXPF(ZZ(7)*C3) 

F SB { C9 ) FMUfZZ+lO). 

STAIZZ+8) 

ZZ(2)=C3*ZZ( 2)*ZZ(7)*ZZ( 1 C ) » ZZ 1 7 ) 

Z Z ( 5)=ZZ(7)*{ 1 . C-ZZ { 8 ) ) 

ZZ( 2 ) =ZZ ( 7 ) — ( ZZ(7)»ZZ(5)-ZZ( 11 ) )/(2.0*ZZ(5)+ZZ{2) ) 
LDA ( ZZ + 2 )’ 

FSB ( ZZ + 7 ) A JP2 ( 19) 

LDA ( ZZ+2 ) 

FMUM. OOOOOl) F SB ( Z Z + 7 ) . 

A JP2 ( 20 ) LDA ( ZZ + 2 ) . 

STA ( ZZ+7 ) SL J ( 15) 

19 Z Z C 2 ) =ZZ ( 7 ) 

BOX 42 

20 LOA(O.O) STAIZZ+5) . 

ST A ( ZZ + 6 1 LDA ( ZZ+2 ) . 

FMU ( ZZ + 2 ) S T A { ZZ + 7) . 

ZZ{ 8)=C2/EXPF{ZZ (2)*C3 )-C9 
LACIZZ+2) F M U ( C 8 ) 

F AD ( 1.0) STA ( ZZ + 4 ) . 

DO 25 1=1,4 

ZZ t 9 ) =ZZ ( 4 ) *CON ( I ) 

ZZ ( 3 ) = ZZ ( 2 ) »ZZ( 2 ) 

ZZ ( 15)=ZZ(3)*ZZ(10)/ZZ(9) 

ZZ( 1 4 ) =Z Z ( 1 1 ) * ( 2 . 0- ZZ ( 9 ) ) 

ZZ C 9 ) = 1 .O-ZZ (9) 

ZZ( 1 2 ) = 1 . 0/ ZZ ( 9 ) 

ZZ(9)=C2/EXPF(ZZ (2)*ZZ(12)*C3)-C9 

Z Z ( 9 ) =W ( I ) /SQRTF { ABSF (ZZ(14)+ZZ(15)*(ZZ(8)-ZZ(9)))) 
FAD ( ZZ + 5 ) ST A ( ZZ+5 ) . 

LDAIZZ+9) FKU(ZZ+12>. 

FMU { Z Z + 1 2 ) FAD ( ZZ + 6 ) . 

STA ( ZZ+6 ) 

25 CONTINUE 

ZZ(7)=2.0*SQRTF( ABSF(ZZ(4 ) ) ) 

ZZ(5)=S*ZZ(5)*ZZ(7) 

ZZ( 6)=ZZ(3)»ZZ(6)*ZZ(7)/ZZ(1 ) 

BOX 43 

ZZ(7)=SINF(ZZ(5) )*C6 
ZZ ( 8 ) =COSF ( ZZ { 5 ) ) 

ZZ(3)=ZZ(5)+ZZ(13) 

ZZ( 4)=SINF(ZZ(3) )*ZZ( 1 J+C.5E-4 
ZZ( 1 0 ) =C 5*ZZ ( 1 )*C0SF(ZZ(3) ) 



BOX 44 



F I ( 1)=XP(2)*XP(6)-XP(3)*XP(5) 
FI(2)=XP(3)*XP(4)-XP( 1 ) * X P { 6 ) 
FI(3)=XP(l)*XP(5)-XP{2)sXP(4) 
ZZ ( 1 )=FI(2)*XP(3)-FI(3)*XP(2) 
ZZ(2 ) = F I ( 3 ) *X P ( 1 > — F I ( 1 ) * XP ( 3 ) 
ZZ( 3)=FI(1 ) *XP ( 2 )-F I ( 2 ) * XP ( 1 ) 
FMU ( ZZ + 3 ) STAUEMP) 

LOA ( ZZ+2 ) FMU(ZZ+2) 

FAD(TEMP) STAUEMP) 

LDA ( ZZ+ 1 ) FMU ( Z Z + 1 ) 

FAD(TEMP) STA(TEKP) 

TEMP=SQRTF( ABSF( TEMP) ) 

L DA ( 1.0) FDV(TEMP) 

STA ( TEMP ) LDA ( 2 Z+ 1 ) 

FMU ( TEMP ) STA ( ZZ+ 1 ) 

LDA ( ZZ+2 ) FMU ( TEMP ) 

STA ( ZZ + 2 ) LDA ( Z Z + 3 ) 

FMU ( TEMP ) STA ( ZZ+3 ) 



BOX 45 



DO 30 1= 1,3 
13=1+3 

ZZ ( 5 ) =XC ( I)*ZZ(8)+ZZ( I ) * ZZ ( 7 ) 

ZZ(9)=B( I , I I )+XC( I > + { B { 13,1 I)+XC( 13) )«ZZ(6) *1.0E14 
DEV ( I, 1 )=ZZ(9)+ZZ(5) 

CEV( I,2) = ZZ(9)-ZZ(5) 

ZZ( 5 ) = XC( I) + ZZ(10)+ZZ(I)«ZZ(4) 

ZZ ( 9 ) =B ( 13, I I ) + XC( 13) 

DEV ( 13, 1 )=ZZ(9)+ZZ(5) 

30 DEV ( 13,2 ) = ZZ ( 9 )-ZZ ( 5 ) 

LDA(C.O),STA(Tl ) ,STA(T2) 

DO 29 1=4,6 

T 1 = T 1 +D6V ( I , 1 )*DEV( 1,1) 

29 T2=T2 + DEV( I ,2 )*DEV{ I, 2) 

DEV (7,1 )=FMASS«T1*0.5180C761500 
DEV(7,2)=FMASS*T2*0.5180C761500 
END 
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BOX 1 



PROGRAM SLAVE 

DIMENSION A8(3) ,e<7,2370),FI(3), I B { 3 ) ,L B ( 237C ) , XVZ ( 3 ) , EMAX ( 30 ) 

Cl MENS ION YENTR V ( 1C ) , ZEN7RY (1 0 ) , CC TA L.(2 370 ) .KET ( 26 ) 

DIMENSION IHIST0(105),NUM8ERS(56).MISTAKE(10) , . 

DIMENSION NSIDE(6),ESIDE(6) ,N I TEMP (10) , IC( 1 00 ) . TSL I ( 30 ) , NTSC ( 3C ) 
COMMON IX, I Y, IZ,FI,A,NLINE,KPLANE,LMAX,IB,LB,XYZ,MAXLAT, B, OCTAL 
COMMON IH ISTO, NUMBERS, AB, I L, I TEMP, ITEM, IT ,TEMP,T£M,KET 
COMMON NSU8ST ,EBUL» ETHER M, Ns IDE ,NBUL , NS A ME, INTERST.INPAI RS ,NFACE 
COMMON ESIDE.NV AC, ECAL,OENERGY,N ITEM P, DEPTH, PERC 1ST, IC,EPB 
COMMON YENTRY , ZENTRY 
EQUIVALENCE (LB, OCTAL) 

REWIND 8 
PRINT 109 

REAO TAPE 8, IX, I Y, I Z , AB ,NL INE .NPLANE , MAXLAT ,LMAX,A , SPHI , FMASS', 

1 ETHRESH, THERMAL, EPB, LULL, ALPHA, BETA, NTS, MNTS . TT IME, TSLI , NTSC, 
20ENERGY,TENERGY, MISTAKE, YENTRY, ZENTRY, NA, EMAX, EFOUNO 

REAO TAPE 8 , ( LB ( L), ( 8 ( I ,L) 1*1,7) L*1,LMAX) 

WRITE OUTPUT TAPE 3,398 

398 FORMAT { 30H SET SLJ SWITCH 2 UP IF NEEDED) 

SLS C 399 ) . 

399 CONTINUE 
PRINT 101 

DO 400 1=1,6 
READ INPUT TAPE 2,100 
100 FORMAT ( 80H 



400 

405 



200 

210 

220 

101 



1 



) 



PRINT 100 

READ INPUT TAPE 2, 405, )( KET ( I ) 1*1,25) 

FCRMATMOX, 1314) 

KET { 26 )= 1 000 
SLJ2( 200 ) SLJ (220) 

PRINT 210 

FORMAT ( 40H SELECTIVE JUMP SWITCH NUM8ER TWO WAS UP) 
PRINT 101 
FORMATt ////) 

IHISTOt 1 )=5H( 13, F 
I H I ST0(2)=4H7. 1 , 

IHI STO( 3 ) =5H3H P. 

IHIST0(4)=6H4HC. 

IHI STO { 105)= IH ) 

I TEMP=32768 

I2E 1 5 = I TEMP 

I2€24=ITEMP«512 

F2E24=I2E24 

I2E30=ITEMP*ITEMP 

MASK25=77770B 

MASK69 =7 777000008 

MASK 13=77770000000008 

MASK14=70000000000000B 

MASKSA=1 OOOOOCCCCOCOOOB 

MASK$U=200000000000000B 

MASK SET =73000000000000078 
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BOX 2 

+ LDA (NTS) AJP1 ( L+2) LDA(MNTS) STA(NTS) 
NSUBST*0 
EBUL=0 • 0 
ETHERM=E FOUND 
DC 1 1105 1*1,6 
ESIDEl I ) =0.0 
1405 NS IDE ( I )*0 
NBUL=0 
NSAME^O 
INTERST*0 
NAA=NA— 1 



DO 1407 1-1, NAA 
1407 TSLI(I)*TSLI(I)*1.0E— 14 
TTfME*TTIME*1.0E-14 
DO 1410 L* 1 » LHAX 
+LDA4 ( LB ) SCL(MASKSET) SlJ2(l+2> 
1410 CONTINUE 



SCL (MASK69 ) SCL (MASKl 3) 



♦STA4ILB) 
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. 





V 



1411 

1412 

1415 

1416 
141? 



1425 



1426 

1427 



1428 

1444 

1445 



1449 



1450 

1451 

1452 



ECX 3 



!:2i 



ICI 

ICI 



( NFA 
(NFA 



1 + 1 

>♦8(7, L) 



DO 1449 L=1,LMAX 
LDQU { L8 ) LDL ( MASK 1 4 ) 

AJP (1415) ALS ( 9 ) 

STA(NFACE) 

NSIOE(NFAC( 

ESI DE (NFACc 
GO TO 1449 
SLJ2( 1416) 

CALL ROUNO(L) 

+ SL J ( L+3 ) 

+LDQ4 ( LB ) LDL (MASK13 ) ARS(27) STA(IL) 

LILKIL) LDA(IL) ALS(27) STA(IL27) 

SLJ 2 ( 1417) LDA4 ( L8 ) SST(IL27) STA4ILB) 
TEMP-BI7.D-THERMAL 
A JP3 ( 1 425 ) RAC(NBUL) 

EBUL=EBUL+B(7.L) 

LCA4 ( LB ) SST(2B) STA4ILB) 
TEM=BJ7»L)-ETHRESH 
A JP 2 ( 1449) 

CALL ROUND ( L ) 

LIL 1 ( IL) 

LDAl(LB) SST ( 4B ) STAKLB) 

GO TO 1449 
ETHERM^ETHERM+BITtL) 

SL J2 ( 1445 ) LOCI (LB ) 

LDL ( 4B ) AJP ( 1 427 ) 

DO 1426 J=1 » L 

LDQ2 ( LB ) LDL( MASK13) 

SUB ( I L27 ) A J P 1 ( 1426) 

LDQ2ILB) LDL ( 28 ) 

AJP1( 1426) 

ITEMP=L 

AL S ( 15) ADD21LB) 

I TEMP=J 

ALS (15) ADD41LB) 

SLJ (1444) 

CONTINUE 
LDAKLB) 

LAC( IL) 

ADD4 ( LB ) 

LDA ( MASKSA) 

LDAKLB) 



SCL(MASKSA) SCL(MASKSU) STA21LB) 
SCL(MASKSA) SCL(MASKSU) ST A4 (LB ) 



SST (4B ) STAKLB) 

INA4 ( C ) AJ P ( 1 428 ) LCA(MASKSU) 

STA4 ( LB ) SLJ ( 1444 ) 

ACD4 ( LB ) STA4 ( LB ) 

SST ( 46 ) STAKLB) SLJ(1449) 



LDAKLB) SST (46 ) STAKLB) 

+LDQ4 ( LB ) LDL ( MASKl 3 ) ARS(27) STA(IL) 

♦ LILKIL) LDAKLB) SST(IB) STAl (L 6 ) 

+ LDL ( MASK69 ) AJPK1449) LAC(IL) INA4(0) 

+ AJPKL + 3) LCA4 ( LB ) SST(NASKSA) ST A4 ( LB ) . SLJ ( 1 449 ) 
+LDA41LB) SST ( MASKSU) STA4 (LB ) 

CONTINUE 

NVAC=0 

DO 1456 L= 1 * LNAX 

LDQ4 (LB ) LDL ( 4B ) AJPK14SC) RAO(NVAC) 

LDQ4 ( LB ) LDL (MASKSA ) AJP(1451) RAO(NSAFE) 

LDQ4(LB) LDL ( MASKSU ) AJP(1452) 

LDQ4 ( LB ) LDL ( MASK69 ) AJP11456) 

NVAC*NVAC+MAXLAT-LMAX 



RAC (NSUBST ) 
RAO (I NTERST) 
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ecx 4 



TEMP=ESI DE l 1 ) 

ITEMP^NSIDE ( 1 ) 

00 1457 1=2,6 
ITEMP=ITEMP+NSIDE( I ) 

1457 TEMP=TEMP+ES IDE ( I ) 

ECAL=ETHERM+EEUL+TEMP 

LMAXCAL= INTERST+NSAME+NSLBST+NBUL+ITEMP 

PRINT 1458,IX*AB ( 1 ) , IY*AB(2) , I Z , AB (3 ) , MAXLAT , A, SPH I , FMASS , ETHRESH, 
1THERMAL 

1458 FORMAT! 29H THE LATTICE DIMENSIONS ARE -5X, 1 1HX-DIRECT ICNI6, 6H 

1 X A =F8. 3, 1CH ANGSTR0MS/34X, 1 1 HY-D IRECT ICNI 6, 6H X A =F8.3, 1 OH 

2 ANGSTR0MS/34X, 11 HZ-DI RECT I ON 16 , 6H X A =F8.3, 10H ANGSTROMS// 21H 

3 THE LATTICE CONTAINS 1 5 , 6H AT0MS//33H THE LATTICE SPACING CONSTANT 

4 A =F6 .3, 1CH ANGSTRCMS/27H THE SPHERE OF INFLUENCE ISF12.3, 1CH 

5 ANGSTROMS/ 19H THE ATOMIC PASS ISF20.3,4H AMU// 

6 40H THRESHOLD ENERGY FCR VACANCY P RCCUCT ICNF8 . 3 , 4H EV. / 

7 25H THERMAL THRESHOLD ENERGYF23. 3, 4H EV.) 

LDA { ALPHA ) FMU(IBO.O) FCV ( 3 . 14 15926 536 ) ST A ( YANGLE ) 

LDA{ BETA ) FMU(180.0) FCV ( 3. 1 41 5926536 ) ST A ( ZANGLE ) 

IT=LMAX-MAXLAT 
PRINT 101 

PRINT 1454, I T, LULL, EPB, ALPHA, YANGLE , BET A, ZANGLE 

1454 FORMAT ( 23H CONDITIONS OF THIS RUN///I3, 20H ATOMS WERE SHOT I Km 

128H ONE ATOM WAS SHOT IN EACHl3,llH TIME STEPS/ 

3 38H ORIGINAL ENERGY CF EACH INCOMING ATCMF12.3,4H EV./// 

4 28H DIRECTION OF INCOMING ATOMS/ 

526H ALPHA = ARC TAN ( VY/VX ) =F8.4,11H RADIANS =F8.3,8H DEGREES/ 
626H BETA = ARCTANI VZ/VX ) =F8.4,11H RADIANS *F8.3,8H OEGRE<€S> 

I TEMP=MAXLAT + 1 

PRINT 1453, ( L ,YENTRY( L-MAXLAT) , ZENTR Y { L-MAXLAT ) L= ITEMP, LM AX ) 

1453 FORMAT ( ///37H POSITIONS CF ENTRY CF INCOMING AT0MS//7H NUMBER,9X, 

1 1HX9X,1HY9X, 1HZ//(I6,8X,5HC.C0C,2F10 .3) ) 

PRINT 1455,NTSfcMNTS,TTIME, t TSLI ( I ) ,NTSC( I ) 1 = 1, NAA) 

1455 FORMAT { 1 Hi , // 41H THE NUMBER OF TIME STEPS OF THIS RUN WASIlO/ 

1 47H THE LIMIT SET FOR THE NUMBER OF TIME STEPS WASI 4/ 31H THE 

2T0TAL TIME OF THIS RUN W AS 1 PE 1 2. 5 ,8H SECONDS// 41H TIME STEP LENGT 
3H NUMBER OF FIRST TIME/4X, 3oHI N SECONDS STEP OF THIS L 

4ENGTH/1/1PE14.4.I18)) 

LDA(NTS) SUB ( 30 ) +AJP2(L>2) LDA ( NT S) STA(IT) SLJI5541 ) 

LDA ( 30 ) STA(IT) 

5541 CONTINUE 
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PRINT 1459,ECAL,TENERGY,CENERGY,ETHERM, ( I,EMAX< I) 1=1, IT) 

1459 FORMAT! 1H1 ,/28H THE ENERGY ACCOUNTED FCR ISF14.3.4H EV./ 

1 28H SUM OF ENERGIES OF ATOMS F14.3.4H EV./ 

2 28H ENERGY SHOT INTO LATTICE F14.3.4H EV./ 

3 45H ENERGY ABSORBED BY LATTICE AS THERMAL ENERGYF1 2 .3, 4H EV.//. 

4 29H NUMBER OF ENERGY OF MOST/ 

5 29H TIME STEP ENERGETIC ATOM// (17 , F 16 .3, 3H EV)) 

I NT ERST= I NTERST/2 

PRING 1460,LMAXCAL,LMAX,NSAME»NSUBST,INTERST,NVAC 

1460 FORMAT! 1H1,/38H THE NUMBER CF ATOMS ACCOUNTED FCR IS 15/ 

1 38H THE TOTAL NUMBER OF ATOMS IS 15// 

2 38H THE NUMBER IN ORIGINIAL POSITION'S IS 15/ 

3 38H THE NUMBER OF REPLACEMENTS IS 15/ 

4 38H THE NUMBER OF INTERSTITIAL PAIRS IS 15/ 

5 38H THE NUMBER OF VACANT LATTICE SITES ISI5////) 

I F( M I STAKE ( 1 ) ) 9955,9955,9956 

9955 PRINT 9957 

9957 FORMAT ! 30H THE ENERGY CHECK DIC NOT FAIL ) 

GO TO 9965 

9956 PRINT 9959, M I STAKE ( 1 ) 

9959 FORMAT! 24H THE ENERGY CHECK FAILED I4,16H DIFFERENT TIMES) 

9965 IF { MI STAKE! 2 ) ) 9961,9961,9962 

9961 PRINT 9963 

9963 FORMAT { 62H NO INTERSTITIALS OR VACANCIES OCCUREO ALONG AN EDGE OR 
1CCRNER > 

GO TO 9966 

9962 PRINT 9960, M I STAKE t 2 ) 

9960 FORMAT ( 60H AN INTERSTITIAL CR VACAICY CCCURRED ALONG AN EDGE OR CO 
• 1RNER I 4 , 1 6H DIFFERENT TIMES) 

9966 CONTINUE 
PRINT 101 
PRINT 1462 

1462 FORMAT! 24 X, 6HNUMBER ,5X',12HT0TAL ENERGY/5H FACE,6X,5HGR0UP,7X, 

18HIN GROUP, 7X,6HIN EV.//) 

PRINT 1464,.NBUL,EBUL, (NSIDE! I ) ,ESIDE 1 1 ) 1 = 1,6) 

1464 FORMAT! 19H KNOCK-CNS I9,F16. 3/19H 1 FRONT FACEIS',F16 

1.3/19H 2 BACK FACE I9,F16.3/19H 3 LEFT FACE I9,F16;3/ 

2 19H 4 RIGHT F ACEI 9 ,F1 6. 3/19H 5 TOP FACE I9,F16.3/ 

3 20H 6 BOTTOM FACE 18 ,F 16.3) 
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BOX 5 



1465 



1466 

1467 
14 r 68 

1469 



1470 

1471 

1472 
1474 



1476 

1480 

480 



PRINT 1 465,NS AME 
FORMAT! lHl»lXtS9H THE NUMBER OF 
SITICNS ISI5//10H THEY ARE/) 

ITEM = 7777 777000077777 B 
IT-77700007777777776 
J — I 1 
K = 1 
L= 1 

DO 1466 1=1.10 

NUMBERS ( I ) = 0 

DC 1480 M=1,LMAX 

IFtJ-11) 1474,1467,1467 

PRINT 1 472*1 NUMBERS ( I ) 1*1,10) 

DO 1468 1=1,10 
NUMBERS! I )=0 
J=1 
L=L + 1 

IFTK-501) 1470,1469,1469 
PRINT 109 
K= 1 
L = 1 

GO TO 1474 

IF(L-ll) 1474,1471,1471 
PRINT 110 
L= 1 

FORMAT! 10110) 

LDA5 ( LB ) SCL(IT) ARS!27) STAUTEMP) 

LAC ( I TEMP) INA510) AJPK 1476) LDA5CLB) 
SCL ( ITEM) AJPK1476) 

NUMBERS! J)=.M 
J=J + 1 
K=K + 1 

IN I 2 (- 1 ) SIL2 ( IT) LDAIIT) AJP1480) 
PRINT 1472, {NUMBERS (I ) 1*1, IT) 

BOX 6 



ATOMS THAT ARE IN THEIR INITIAL PO 



PRINT 1560, NVAC 

1560 FORMATdHl, 1X,37HTHE NUMBER OF VACANT LATTICE SITES ISI6///10H THE 
1Y ARE //) 

IT-1 

DO 1565 L=1 , MAXLAT 
L0Q4 ( LB ) LDL14B) AJPK 1565) 

NUMBERS! IT )=L 

RAOIIT) SUB (51) AJPK1565) 

PRINT 1472, (NUMBERS! IT) IT-1,50) 

IT-1 

CONTINUE 

PRINT 1472,(NUMBERS(I ) I«1»IT) 



1562 



1565 

1566 
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BOX 7 



1499 

1500 

1502 



1515 

1517 

1518 



1520 



PRINT 109 
PRINT 1499 

FORMAT ( 1 3H REPLACEMENTS//) 

PRINT 1 500,NSUBST 

FORMAT! IX, 30HTHE NUMBER OF REPLACEMENTS 
1 INITIAL FINAL FENETRAT ION 

DO 1502 1=1,100 • 

IC< I ) =0 

DO 1517 L= 1 » LMA X 

LDQ4ILB) LDL ( MASKSU) AJP(1517) LDL IMASKl 3 ) 
CALL INITIAL(L) 

CALL OEPTHP(L) 



IS 16/// 
PENETRATION IN/ 



ARS ( 27 ) STA(IL) 



PRINT 1 51 5,L, IL» DEPTH. PERDIST 
FORMAT! 17,1 10, FI 6. 3, FI 8. 3) 



1521 

1525 

1527 



1530 



CONTINUE 
PRINT 109 
PRINT 1499 
PRINT 1520 

FORMAT 1 40H DISTRIBUTION CF X-OIRECTICN PENETRATICN//2 IX, 

117HR A N G E18X,6HNUMBER/) 

DO 1521 1=1,30 
J = I-7 
TEM- J 

TEM=TEM-0.5 
TEMP = TEM+;i.O 

PRINT 1525, J. TEM, TEMP, IC(IJ 

FORMAT! 13, F11.1.9H X A - 16HX PENETRATION -F6.1,4H X AI11) 
DO 1527 1=1,30 
NUMBERS! I )=>IC!I) 

IL=-7 

CALL HISTOI30) 

PRINT 109 
PRINT 1499 
PRINT 1530 
FORMAT l 35H 
117HR A 



DO 

IT' 



1531 

1-1 



DISTRIBUTION CF RADIAL PE N6TR ATI 0N//2 1 X, 
N _ G E22X.6HNUMBER/) 



1 = 1,20 



1531 PRINT 1532,1 , IT, I, IC< 1+30) 

1532 FORMAT! 13, 19, 9H X A - 21HRA0I AL PENETRATION -14, 4H X AI12) 
DO 1533 1=1,20 

1533 NUMBERS! I )=IC! 1 + 30) 

CALL H ISTO { 20 ) 



S5H 

56H 
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BOX 8 



PRINT 109 
PRINT 1534 

1534 FORMAT ( 1 9H INTERSTITIAL PAIRS//) 

PRINT 1535,INTERST 

1535 FORMAT! 10H THERE AREI4.19H INTERSTITIAL 
1 THE PAIR NUMBER PENETRATION OF 
2PENETRAT ION OF LL 



3 L LL OF 

4 IN X-DIRECT ION 
DO 1538 1*1,100 

1538 ‘ICC I ) =0 

DO 1550 L=; 1 , LMAX 
LDQ4 ( LB ) LDL (2B ) 



PA IRS// 1 1H NUMBERS 0F/TC4H 
L PENETRATION OF L 

PENETRATION OF LL/ 1GSH 

SITE. IN X-DIRECTICN IN RADIAL OIRECTION 
IN RADIAL OIRECTION//) 



AJP 1 ( 1550) 



LDQ41LB) LDL ( MASK69) AJP (1550) ARS115) 

ST A ( LL) LDL ( MASK 13 ) ARS!27) STA(IL) 
LILKLL) LDQML8) L0L12B) AJPK155C) 
CALL INITIAL ( L) 

CALL DEPTHP(L) 

TE =DEPTH 
TEM = PERD I ST 
CALL INITIAL(LL) 

CALL DEPTHP ( LL ) 

PRINT 1545,L, LL, IL . TE , TEM , OEPTH, PERD I ST 
1545 FORMAT! 15,16, 19, F16.4,F2C.4,F22.4,F2 1 .4 ) 
1550 CONTINUE 

DO 1555 1=1, ICO 

1555 ICC I) = IC(I )/2 
PRINT 109 
PRINT 1534 
PRINT 1520 

DO 1556 1=1,30 
J = I - 7 
TEM = J 

TEM=TEM— 0.5 
TEMP=TEM+:i.O 

1556 PRINT 1525, J, TEM, TEMP, IC (I ) 

CO 1557 1=1,30 

1557 NUMBERS! I ) = IC( I ) 

IL=-7 

CALL HIST0130) 

PRINT 109 
PRINT 1534 
PRINT 1530 
DO 1558 1=1,20 
I T= I— 1 

1558 PRINT 1 532, I,IT,I,IC( I+3C ) 

DO 1559 1=1,20 

1559 NUMBERS! I)=IC! 1+30) 

CALL HISTOI20) 
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BOX 9 



1571 

1572 

573 

574 

1576 

1577 



L7X,1 HX8X, 1HY8X, 1HZ8X , 2HVX7X, 2HVY7X 
11HX RADIAL//) 



1567 



1578 



1569 

1579 



1573 

1574 



1580 



1680 

1526 

1681 



LDA(NBUL) A JPl 9951 ) 

PRINT 109 
PRINT 1572 

FORMATl 10H KNOCK-ONS//) 

PRINT 573 » N8UL 

FORMATl 10H THERE ARE I5t 1 CH KNOCK-ONS//) 

PRINT 574 

FORMAT 191X»1 1HPENETRATI0N/4H 
1 ,2HVZ5X,6HENERGY7X,5HLB( D1CX, 

FORMATl I5,7F9.3,3X,015,2F8.2) 

DO 1577 1=1,100 

ICCI )=0 

IT»1 

00 1578 L=1 » LMAX 
LDQ41LB) LDL128) A J P t 1578) 

LDA(IT) SUB ( 5 1 ) AJPH 1 567) 

PRINT 109 

PRINT 1572 
PRINT 574 

1 T= 1 

CALL INITIAL (L) 

CALL DEPTHP(L) 

CALL ECOUNT(L) 

PRINT 1576.L, (B( I,L) 1=1 ,7 ) .OCTAL l L) .DEPTH, PERDI ST 
RAOl IT) 

CONTINUE 
PRINT 109 
PRINT 1572 
PRINT 1520 
DO 1569 1=1,30 
J = I -7 
TEM= J 

TEM=TEM-0. 5 
TEMP=TEM+1 .0 

PRINT 1525 t .J,TEM,TEMP,IC (I L 
DO 1579 1=1,30 
NUMBERS ( I )=IC ( I ) 

I L*-7 

CALL HISTOl 30 ) 

PRINT 109 
PRINT 1572 
PRINT 1530 
DO 1574 1=1,20 
I T= I— 1 

NUMBERSl I ) = IC ( 1+30) 

PRINT 1532..I, IT, I, IC l I+3C) 

IL“0 

CALL HI STO (20 ) 

PRINT 109 
PRINT 1572 
PRINT 1580 

FORMAT (23H DISTRIBUTION CF ENERGY// 

115X17HR A N G El 4X, 6HNUM8ER/ ) 

DO 1680 1=1,25 

PRINT 1 526 , I , KET ( I ),KET{ 1+1) ,IC( 1+50 ) 

FORMATl 13,17, 21H EV - ENERGY -16, 3H EVI10) 
CO 1681 1=1,25 
NUMBERSl I )=uC ( 1 + 50) 

CALL HISTOl 25 ) 
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ecx 10 



9951 DO 1599 J = 1 ,6 

IF IJ-l) 1581,1581,1583 

1581 PRINT 1 582»NS IDE!1) 

1582 FORMAT! 1H1 , IX , 1 6HSPUTTERED 
GO TO 1585 

1583 PRINT 1584,J,NSI DEI J) 

1584 FORMAT! 1H1, IX, 45HTHE NUMBER 
1H ISIS) 



ATOMS-I 10 , 22H ATOMS WERE SPUTTEREO.) 



OF ATOMS WHICH EXITEO THROUGH SIDE-12, 3 



1585 LDA21NSIDE) AJP(1599) AJP311599) 

PRINT 1586 

1586 FORMAT! ///91X, 23HPENETRAT ION TIME STEP/4H L7X, 1HX8X, 1HY8X , 1 
1HZ8X , 2HVX7X ,2HVY7X,2HVZ5 X,6HENERGY7X,5HL8(L) 10X,22HX RADIAL 
20F EXIT//) 

1587 FORMAT! 15 ,7F9 .3 , 3X, 01 5 ,2F8.2 , 19) 

NFACE=0 



1 = 1,100 



L= 1 



1594 



1695 



, UMAX 

LDHMASK14) 
S IL2 ! IT) 

ARS ( 3 ) 



DO 1589 
1589 (1C! I ) =0 
NTS = 0 
DO 1595 
LDQ4 ( LB ) 

STA(NFACE) 

LDL IMASK25) 

CALL INITIAL !L) 

CALL DEPTHP(L) 

CALL ECOUNT(L) 

PRINT 1587, L, IB! I, L) 

1595 CONTINUE 

PRINT 1695, J 
FORMAT! 1H1,12H EXITED 
PRINT 1520 
DC 1696 1=1, 3C 
L = I -7 
TEM = L 

TEM=TEM-0. 5 
TEMP=TEM+ 1 • 0 
NUMBERS! I)=IC( I ) 

PRINT 1525, L,TEM, TEMP, 1C II) 
I L=-7 

CALL HI STO! 30 ) 

PRINT 1695, J 
PRINT 1530 
DO 1697 1 = 1 , 20 
I T= I— 1 

NUMBERS! I ) =, I C C 1 + 30) 

PRINT 1532, I, IT, I , I C ( I +30 ) 
IL-0 



A JP ( 1 595 ) 
SUB! IT) 
STA(NTS) 



ALS ( 9 ) 
AJP1( 1595) 



1 = 1 ,7), OCTAL! L) ,OEPTH, PERO I ST, NTS 



FACE 14 // ) 



1596 

1696 



1597 

1697 



CALL H I STO ( 20 ) 

PRINT 1695, J 
PRINT 1580 
DO 1698 1=1 ,25 

1598 NUMBERS!I)=IC(I+50) 

1698 PRINT 1526,I,KET!I),KET! I + 1),IC!I + 50) 
CALL HIST0125) 

1599 CONTINUE 
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ecx li 



107 

108 
109 

m 
112 
1 13 
1 14 

115 

116 

117 

118 

119 

1 

120 



SL J 1(120) 

PRINT 119 
L = 1 

PRINT 118.L, (8< I,L> I** 1 , 7 ) .OCTAL ( L ) 
FORMAT ( 1H1) 

FORMAT (/) 

DO 117 Li2»LMAX 

IFC 1 0 -» ( { L — 1 1/101-L+l ) 117,113,113 
IF( 50*( (L-l )/5C)-L+l) 116,114,114 
PRINT 119 
GO TO 117 
PRINT 110 

PRINT 118,L,(E(I ,L) 1=1 ,7),CCTAL(L) 
FORMAT ( IS.7F9.3,3X,015) 

FORMAT (1H1»///13H LATTICE DUMP///4H 
,2HVY7X,2HVZ5X,6HENERGY 7X,5HL8(L)/J 
ENO 



L7X, 1HX8X, 1HY8X, 1HZ8X,2HVX7X 



BOX 12 

SUBROUTINE ROUNO(L) 

CIMENSION A B ( 3 ) »B(7 ,2370 ),FI (3) , 18 (3 ) , LB ( 2370 ) , XYZ ( 3 ) , EMAX ( 30 ) 
DIMENSION YENTRY ( 10), Z ENTRY! 10 ) ,OCTAL( 2370) ,KET(26) 

DIMENSION IHIST0C105),NUMBERS(50),MI STAKE (10) 

DIMENSION NSIDE16) ,ESIDE (6) ,NI TEMPI! 0) ,IC(100),TSLI (30),NTSC( 20) 
CCMMCN IX, I Y, IZ,FI , A ,N L I NE , NPLANE , LM AX , I 8 ,L 8 ,XYZ , MAX LAT, 8, OCT AL 
COMMON I HISTO, NUMBERS, A8,I L, I TEMP, ITEM, IT, TEMP, T EM, KET 
COMMON NSU8ST , E8UL, ETHERM ,NS I DE , N8UL ,NSAME, INTERST, INPAIRS ,N FACE 
COMMON ES IDE, NV AC, EC AL,OENERGY,N ITEM P, DEPTH, PERCI ST, IC.EP8 
COMMON YENTRY ,ZENTRY 
EQUIVALENCE (LB, OCTAL) 

DO 1320 1=1.3 
1320 IBI I )=B( I,L)/A+C.5 

ITEMP=IB( 1 )+I8(2)+I8(3) 

IF.(2*( ITEMP/2J-ITEMP) 1 330 , 1325, 1 330 
1325 IL=IB(3)*NPLANE+I8(2)*NL INE+I8! 1 )/2+ 1 
RETURN 

1330 NUMBERS ( 1 ) ^ 18(3) *NPL ANE+ 18(2) *NL INE + ( IB ( 1 ) + 1 ) /2+ 1 

NUMBERS (2 )= I E ( 3 ) *NPL ANE ♦ 18(2) *NL INE+ ( I B ( 1 )-l ) /2+ 1 

NUMBERS ( 3 )= 18(3) *NPL ANE + ( I 8 ( 2 ) + l ) «NL INE+ 18(1) /2+1 

NUMBERS ( 4 ) = IB(3) *NPLANE+ ( I B ( 2 )— 1 )*NLINE+ IB ( 1 ) /2+1 

NUMBERS! 5)=(I8(3)+1 )*NPLANE+ 18(2) »NL INE+ IBM) /2+1 

NUMBERS(6)=(I8(3)-1 )*NPLANE+ 18(2) *NL INE+ 18(1) /2+1 

TEMP= 1 .0E50 
DO 1340 1=1,6 
IT-NUMBERS! I ) 

AJP311340) SUB ( M AXL AT ) 

A JP { 1333) AJP2( 1340) 

1333 CALL INITIAL(IT) 

TEM=0. 0 

DO 1335 ITEMP^l , 3 

Fit ITEMP )=XYZ( ITEMP )— BC I TEMP,L ) 

1335 TEM=TEM+FI( ITEMP) *FI ( ITEMP) 

IF! TEM-TEMP ) 1337,1340,1340 
1337 IL=IT 

TEMP=TEM 
1340 CONTINUE 
END 
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BCX 13 



SUBROUTINE INITIAL (L) 

DIMENSION AB(3) , 8 ( 7 , 2 370 ) , F I (3), I 8(3) , L8 (237C ) ,XYZ ( 3 ) 
DIMENSION YENTRYI 10) , ZEN TRY! 1 0 ) , CCTA L ( 23 70 ) , KET ( 26 ) 
DIMENSION IHIST0I105) 

NSI0EI6) 



EMAX ( 3C ) 



DIMENSION 
COMMON 
COMMON 
COMMON 
COMMON 
COMMON 



05) , NUMBER SI 50) .MI STAKE! 10) 
, ESIOE (6 ) , NITEMP (10) ,IC(10C 



) ,TSLI( 30) , NTSCI3C) 



IX, I Y, I Z, FI » A ,NLI NE»N'PLANE» LM AX , I B, LB, XYZ » MAX LA T, 8, OCTAL 
I H IS TO, NUMBERS, AB ,IL,ITEMP,ITEM,IT, TEMP, T EM ,KET 
NSUBST ,EBUL, ETHER M , NS IDE , NBUL , NSAME, I NT ERST, INPAIRS, NF ACE 
ESIDE,NVAC,ECAL,OENER GY, NIT EM P, DEPTH, PER CIST, IC,EPB 
YENTRY , ZENTRY 
EQU I VALENCE (LB,OCTAL) 

IFIL-MAXLAT) 1003,1003,1023 
1003 CONTINUE 

IBt3)=(L— 1)/NPLANE 
STQ ( ITEMP ). 

IB(2)=ITEMP/NLINE 
STQ( ITEMP). 

1 8 { 1 )=ITEMP«2 
ITEMP=IB(2)+IB(3) 

IFC 2* ( I TEMP /2)- 1 TEMP) 1 0C5 , 1 0 1 0 , 1 0 10 
1005 'I B 1 1 )=IB( 1 ) + 1 
1010 DO 1 C 1 5 1=1,3 
F I ( I ) = I B ( I ) 

1015 XYZ ( I ) =A*FI ( l ) 

RETURN 

1023 XYZ ( 1 )=0.0 

XYZ (2)=YENTRY(L-MAXLAT) 

XYZ ( 3 ) = ZENTRY ( L-MAXLAT ) 

END 



BOX 14 

SUBROUTINE ECCUNTIL) 

DIMENSION AB ( 3 ) , B ( 7 .2 370 ) , F I (3), 1 8(3 ) , LB ( 2370 ) , XYZ ( 3 ) , EM AX ( 30 ) 
DIMENSION YENTRY( 10), Z ENTRY! 1 0 ) , OCTA L ( 2370 ) , KET (26) 

Cl MENS ION IH I STO ( 105) , NUMBER S ( 50 ) , MI STAKE ( 10) 

DIMENSION NS ICE(6) ,ESICE(6) ,NITEMP(1 C) , I C ( 1 00 ) , TSL I ( 30 ) , NTSC ( 30) 
COMMON IX,IY,IZ,FI, A, NL I NE , N PL AN E, LM AX , I B ,LB , XYZ , M AXL AT, B , OCT AL 
COMMON IHISTO, NUMBERS, AB,IL, I T EMP, IT EM , I T , TEMP , T EM ,K ET 
COMMON NSUBST,EBUL,ETHERM,NSIDE,NBUL , NSAME, I NTERST.INPA IRS * NFACE 
COMMON ESIDE, NVAC,ECAL»OENERGY, NITEMP, DEPTH, PEROIST, IC,EPB 
COMMON YENTRY, ZENTRY 
EQUIVALENCE (LB, OCTAL) 

TEMP=KET (25 ) 

TEMP=B(7»L )— T EMP 
EN I 1(25) AJP2(25) 

DO 20 1=1,24 
TEMP=KET ( 1+ 1 ) 

TEMP=TEMP-B(7,L) 

A JP2 ( 25 ) 

20 CONTINUE 
25 RA0KIC + 50) 

END 
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ecx 15 



1625 

1630 

1631 

1635 

1640 

1645 

1650 



DIMENSION 
DIMENSION 
DIMENSION 
COMMON IX, 
COMMON 
COMMON 
CCMMCN 
COMMON 



SUBROUTINE DEPTHP(L) 

DIMENSION A B ( 3 ) »B(7 »2370) , F I (3 ) » I B ( 3 ) » L8 (2370 , XYZ ( 3 ) , EM A X ( 3 G > 
YENTRY ( 10) «Z ENTRY! 1 0 ) , OCTA L ! 23 70 ) , KET { 26 ) 

IH I STO (105), NUMBERS (50), MI STAKE! 10) 

NSIDE (6) ,ESI DE ( 6 ) » N I TEMP ( 1 0 ) » IC( 100 ) ,TSL I (30 ) , NTSC! 30 ) 
r IY , IZ.FI ,A,NLINE,NPLANE,LMAX, I B , L B , XYZ , M AX LAT, B, OCTAL 
IH I S TO, NUMBERS »AB ,IL, ITEMP, IT EM, IT , TEMP.T EM , KET 
NSUBST »EBUL»ETHERM,N$IDE»NBUL ,NSAME, INTERST, INPAI RS.NFACE 
ESIDE,NVAC,ECAL,OENERGY,NITEMP,DEPTH,PERDIST, IC,EPB 

YENTRY, ZENTRY 

EQU I VALENCE ( LB »OCTAL ) 

DEP TH=B ! 1 » L )-XYZ ! 1 ) 

1 = 1 

TEM P= I 

IF! (TEMP-6. S)«A-DEPTH) 1630, 1635,1635 
IF.( 1-29 ) 1631,1631,1635 
1 = 1 + 1 

GO TO 1625 
iIC< I ).= IC ( I ) + 1 

PER0IST=SQRTF<(B(2,L)-XYZ(2) )**2+(B( 3,L)-XYZ(3) )**2) 

1 = 1 

T EMP= I 

IFCPERDIST-TEMP + A) 1650, 1645,1645 
1 = 1 + 1 

GO TC 1640 

IC ( I +30 ) = IC ( I *30 > + 1 

END 



EOX 16 

■SUBROUTINE HI STO (NGROUPS ) 

DIMENSION AB(3) , 8(7,2370 ),FI(3),I8(3), LB (2370) , XYZ(3) , EM AX ( 30 
DIMENSION YENTRY (10), ZENTRY ( 1 0 ), OCTAL ( 2370 ) ,KET (26 ) 

DIMENSION I HI STO ( 105) ♦ NUMBERS (50 ), MI STAKE! 10) 

DIMENSION NS I DE ( 6 ) , ES I DE (6 ) , NITEMP( 1 C ) , I C ( 1 CO ) , TSL I ( 30 ) , NTSC ( 3C ) 
CCMMCN IX, I Y, I Z, FI ,A , NL I NE , N PL ANE ,LM A X , I B, LB ,XYZ , MAX L AT, E , OCT AL 
COMMON IHIST0,NUM8ERS„AB, I L, ITEMP, ITEM, I T ,T EMP , TEM, KET 
COMMON NSUBST »E8UL»ETHERM,NSI DE »N8UL , NS A ME, INTERST, INPAI RS ,NFACE 
COMMON ES IDE, NV AC, ECAL,0 ENERGY, NIT EMP, DEPTH, PERCI ST, IC,EP8 
• CCMMCN YENTRY, ZENTRY 
EQUIVALENCE! LB', OCTAL) 

KT0TAL=0 
NUMB I G=0 

DO 1110 1=1, NGROUPS 
NTOTAL=iNTOTAL+NUMBERS( I ) 

IF( NUMBERS ( I )-NUMBIG) 1 1 10 , 1 1 1 0 , 1 1 08 

1108 NUMB I G=NUMBERS ( I ) 

1110 CONTINUE 

PRINT 1109, NTOT AL 

1109 FORM AT ( / / 1 0X , 1 7H 1 00 PERCENT ME ANSI 5/ ) 

TCTALN=NTOT AL 

B I GNUM=NUMB I G 
DO 1113 J=1, NGROUPS 
FNUMBER=N UMBERS ( J) 

NUMBER= 1 00 . C+FNUMBER/B IGNUM+O .5 
TEMP = 1 00. 0«-FNUMBER / TOTAL N 
DO 1111 1=5,104 

1111 ‘I H I S T C ( I ) = 3H1 H 
ITEMP=NUMBER+4 

DO 1112 1=5, ITEMP 

1112 IH ISTO ( I ) =3H 1HX 
IT= J+I L 

1113 PRINT IHISTO.IT, TEMP 
END 

END 



97 



PROGRAM RON 

DIMENSION AB(3) ,8(3,2370 ),NT (2000 ),P (3) ,X!3,2000) ,S 
DIMENSION 18:3) , FI 13) ,XYZ (3 ), NUMBERS ( 6 ) , YENTRY ( 10), 
+ SLJKL + 1) SL J ( 4 ) 

IT**— 1 

WRITE TAPE 6, IT, IT, IT, IT 
4 REWIND 6 
PRINT 3 
3 FORMAT! 1H1 ) 

PRINT 1 

1 FORMAT!////) 

DO 5 1=1,6 
READ INPUT TAPE 2,2 
‘ FORMAT ( 80H 



S( 3) 

ZENTRY! 10) 



1 



) 



5 PRINT 2 
PRINT 1 
SLJ2I7) 

RE AO TAPE 6, IX, IY, IZ, AB, MAXLAT , A , SPH I , FM ASS , ETHRESH , THER MA L , MNP 
READ TAPE 6,((B(I t J) 1=1,3) J= 1 ,MAXL AT ) , NLI NE, NPLANE 
PRINT 1098, IX, AB! 1 ) , I Y ,A6 ( 2 ) , I Z , AB ( 3 ) , MAXLAT , A,SPHI , FMASS, ETHRESH, 
1 THERMAL 

1098 FORMAT { 29H THE LATTICE DIMENSIONS ARE -5X , 1 1 HX-D IRECT I ON I 6 , 6H 



2 

3 

4 

5 

6 
7 



_ _ 1 UIN A U , 

X A = F 1 1.7,1 OH ANGSTR0MS/34X, 1 1HY- C I RECT I ON 1 6 . 6H X A =F11.7, 1 CH 
ANGSTROMS/34 X, 1 1HZ-D IRECT I CNI6 » 6H X A =F11.7,l0H ANGSTROMS// 21H 
THE LATTICE CONTAN I S 1 5, 6H ATCMS//33H THE LATTICE SPACING CONSTANT 
A =F 1 0 • 7, 1CH ANGSTR0MS/27H THE SPHERE OF INFLUENCE ISF16.7, 10H 

ANGSTROMS/ 1 9H THE ATOMIC MASS ISF24.7,4H AMU// 

4 CH THRESHOLD ENERGY FOR VACANCY P RCDLCT I CNF 1 0 . 3 ,4H EV./ 

25H THERMAL THRESHOLD ENERGYF25. 3, 4H EV.) 

PRINT 1 
LMAX=MAXLAT+MNP 

7 READ INPUT TAPE 2,8,T 

8 FORMAT(20X,F10.6) 

DX=T»A 

RA= 1 . 0/A 
ENI6( IB) 

10 READ TAPE 6,IT,P 
LDA(IT) AJP3! 30) 

CO 15 1=1,3 
T = ABSF(P< I)-B('I, IL)) 

FSB(DX) AJP2 (20 
15 CONTINUE 
SLJ (10) 

20 DO 25 1=1*3 
B( I , IL)=P( I ) 

25 X ( I ,N)=P( I) 

LDA(IT) ST A6 ( NT ) +ISK6(3000) SLJ(IC) ENI6(3001 ) 

30 'IN 1 6 (- 1 ) , S I L6( I TEMP ) ENI2(0B) 

DO 45 L=1 ,LMAX 

....... ST A ( NO) 



LIU5 ( IT ) SIL5UN) LIL4UT) SIL4(IL) 



ENA(-rl ) 



LIL1 ( IU) 
SIL1 (IN) 



SI L 1 ( I I ) 
LCA(NO) 



LDA ( 1 1) 
AJP3! 32) 



SUB ( IL ) 

♦ I SK5 ( 4 ) SL J ( 4 2 ) 



E N 1 5 ( 0 B ) S I L4 ( I L ) 

CO 4 5 K= 1,1 TEMP 
LDA3 ( NT ) STA(IU) 

AJPK45) LIUKIU) 

PRINT 41 

41 FORMAT!/) 

LDA ( NO ) A JP2 ( 42 ) 

32 PRINT 35, IL 

35 FORMAT ( // IX ,20H0RBIT OF ATOM NUMB ERI 5// 

1 1 OH NUMBER 0F9X , 1 5HACTUA L PC S I T I ON 16 X , 9HNUMB ER 0F7X , 1 9HNCRMAL I ZED 
3POSITICN/10H TIME STEP6X , 1 HX9X , 1 HY9X , 1 HZ 1 3X , 9HT I ME STEP6X, 1 HX9X, 
4lHY9X,lHZ/) 

IN 1 2 ( 6) 

N0= 1 

42 DO 40 1=1,3 

40 S ( I )=RA*X( I , K ) 

PRINT 43.IN, ( Xt I ,K) 1=1 , 3 ) , IN , ( S ( I ) 1=1,3) 

43 FORMAT! 1 7 , F 1 3 . 4 , 2F 1 0 . 4 , I 16 , F 1 3. 4 , 2F1 C .4 , 19 ) 

45 CONTINUE 

LDA ( IT ) AJP3(5C) ENI6(1B) SLJ(IO) 

50 CONTINUE 
END 
END 
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APPENDIX IV 
PROGRAM FLOW CHARTS 




SYMBOLS USED IN THIS APPENDIX 



START 



a program ^ subroutine 





Read in or write out 




Decision or branch 



General other computer operation 




Position (designation 



[uk-lmnI 

(A) 

(Q) 



132 



Logical product ITK and LMN 
Contents o*f the " A register ss 
Contents erf the n Q. register 

Statement number 7S2. 



, ^ 
age 
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PROGRAM MASTER 



Boxes I B 




100 



Boxes (o-\ 0 




101 




* e>ox 12. 



102 




Boxes \Z-!5 




103 









Boxes 14-17 





104 




Bo* 17 




DO MS 
K r I ,KMAX 



ITEfAP= NM6<) 




105 









I 



Box 18 








106 




Bo x 1 9 




107 











Box 20 




TEMP* % cold - tY 1 " 



i 



TEMprJ 

&(4’W 


12 


. 5 IS ■ • • ^MASS 'TEW? 
l)=TeKP‘C0LD(4-fe) 
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! 



box 21 




IN3=0 



Go To 277 
(Box -25) 



TLAST(M) * TTIME-TSL+TIME 
42 2(C)- 10 1 * PERCENT 
TEMP = B(7> N\) 

LOCATE =0 




Go To 513 

(another page) 



Clear bit *2 of LbCM> 
ZTEMP = lO-ri octal positions 
TTtMP=6-4 octal positions 




Clear MASK2 
bits of LB(M) 



TT£MP= 




Set bit 4*2 


TTEMP-Sf' 5 




of LBCM) 



Clear bit **3 
of CBtlTEMP) 



Clear MASKS bits 
of- L6(M) 

ITEM? =6 it 43 of 
UZC^ • C\eer 
bit 4S of UCM) 



Clear MASK I blfc 
of LfctTTEMP) 



Go To 514- 

(anoihff page) 



CATER.(W2)= M 
W2.=WH+ I 



MISTAKEN) = 
MrSTAKEC3)4 1 







& 









Clear MASK! 
bits of 
LBCM) 





GOO 



l 



MJ - M 
MOVE = I 
NOW =0 
CALL Tk/VAC (mt) 
MOVE -O 
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Go To 


(a) 


Add one to 


k)e*t 




MISTAKE (5) 


?a^e 








J 



Go To 5 \ 1 > 
(next pa<je) 



Put TTEMP in 


* 


* 01 


lb(m) 




ITEMP=I-2' 5 

JTEMP = M-3. ,? 




I = NA MV 


Put JTEMP in 






LB(NAV)^)-ELBCC> 







N ext 



pacje 
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Ill 





e>o* 22 
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Boxes ,23-25 
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Boxes 2C>~2S 




oo *05 
J=|,IMMAX 






<o 




ENERGY = ZCl.K) 


'ENERGY- >-* 


IEWEM.Y* K ■ 








i. 


>o 




305 









BCl-3,K>«B(l-3,K^ 
+3(4-4, K> TSL 


_ \ 


[ 


Set bit *1 

MASK4 oP 


and cleat* 

LB(K) 



N\T= LATERClT) CALL INVACCMJ) 
Upper hal-f c>4 IU * M 
Loner hal-f & TO = M J 
Write Tape 6 , XU > 6 



>o 







DO 241 








ITI =KI2“ l 






261 




Go To 




ir= uiti 






Boy 











t>0 79^ 8 




DO BOG 


T=1,INWAX 




K=I,XT| 




Put KJFACE m octal *1+ of LVjLtti) 




799? 



NFACt= 0 



ITEMP= octal 10-13 of LB(fA') 

TTEhP= octal C-9 of LB(M') 
Clear MA$<* bt U of Lfc(M> 



Clear 4B bit 




C\ear MASK\ blH 


oP LBdTEMP'i 


""XTTEMP >- — ► 


oP LB(JTEMP) 



*0 
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Box 30 
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SUBROUTINE SIT£(M) 



Box 31 




117 






SUBROUTINE INVAUMJ) 



2>o x 32 
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Box 33 




DO 2115 


y«\. 


>0 


JA=TAH 


1 = 1, 3 


►<IPACE(I)> 




TW=TA“l 



£0 



=o 



G> 




*0 


JC=I 










1 


k 




MIST A KE( 2 \ - MI VTAKE( a - ) +- 1 


“ 





TPMAX= S’ 




<2 




TPM AX = ^ 


JPMAX=6 


1 






r \ 




1 









r 


JPMAX=I3 



N00M(moVo 
H0ONi(l-\S)=O 
1= 1 
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(A) 


wl l\ cowtai n 


2.451 


IL4-NLIME + TL 


2452 


IL-NLlbiE + TL 


2453 


XL + NPLANE + TL 


2.454 


XL- NPLANJE 40-L„ 


2455 


IL p 


2456 


XL- WLIKJE4TL+1. 


2457 


XL4KJPLAME4JL-H 


245S 


XL+ NPLAKI.E4NLXNE 


245<* 


JL4WPLANE-NLXNE 


2460 


IL -NPLANE+TL+I 


2461 


1L-NPLANE4NLXNE 


2462 


il-nplane-nlin/e 


2465 


XL~\ . . • 


2470 


IL4NLINE 43L + I 




n\ookj(x)= itemp 
r - i + i 



TPMAX = I- 


-1 














Jp=l 



Noonctp)= 

MOONOO 



=o 







DO 2485 
K=I,TPMAX 




3 


r 




X = 

MOOWCkA 
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Box 3^ 
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subroutine inter. (mq,ii,s) 



Boxes 40-41 
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Box- -4E> 
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PROGRAM SLAVE 



Box I 




START 




REWIND TAPE 8 



3 



/ f>£/\D TAPE 8 

ix, iy,i2, abo-3),n/lia/£, mplame.,maxlat,lmax, 

A, SPHI, FMASS, ET/-/R£SH» THERMAL, EP& 0 LULL, 
ALPHA.* SETA, UTS, MN/TS , TTIME ,TS Li Cl -30) , 
A/TSCO-So'), OBKJERGY , TEV£R6,y , MISTAKE n - 10), 
YEA/TRY Cl- 10} , ZEWTRY 0-103 , UA , EMAX( l- 30>, 
\^EFOb>MD > LBC l-LMAY), BQ-7, l-LW) 



Tupe out veminciev to 
set selective switch 2. up 




a 



RE AO 1 CARO (TAPE 2 ) 



0* 



DO 4-oo 

r* »,& 



PRINT i LINE 
(COMMEKJT READ IN) 



■Roo 




READ XSY KETC l-^s) 
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130 




Box 2 
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Sox 3 




* 



0*i 



av\ot>hev 



pa^e 
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* Next pa^e 
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1445 



I 



Set bit *=4 of Lg Cl) 






Go "To 


^ *> s' 










1 






Set MASKS 0 






bits of 


i-BO-} 





► x= 


IL 


1 


r 


Set bit*bl 
of LB(I) 




















&©x 

4 




NVAC=NVAC+MAXUT- LMAX 




1456 




Add one to 

IMTERST 
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Box 4". 




( PRINT: IX, ASO), 

V SPMI, , ETMR 
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Box G 
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Box 7 




= 0 




1 r 

-JL 



1517 



1521 



PRINT: T^TEM, 
TEMP, IC(l) 





i 



WUMBERS( l-30^=IC(l-30) 




I L* ~7 





CALL HISTOfeo) 



1531 



PRINT: I, IT, 
I, TC(X-*- 36 s ) 




DO 


ISAI 


1 = 


1 > 30 



7= 1-7 
TEM = 7-0.5 
TEMP* 7+0.5 





DO 1531 

x* i , ao 




1 






IT =1-1 
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Box 8 
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»{VRINT: NEUL? 



7 - 1-1 
TEM- J-O.S 

TEfAP = J+O.S 



f PRINT: LjE(l-7 ) iX > 
I OCTALO-^ DEPTH, 
V PERDTST 



IC(HOC»>=0 
IT = l 




Box 9 



DO \5-7? L=t)LMAX 



Skip to next \ 

pa<je of ouipuft J 
papev' J 




CALL INITIALS) 
CALL DEPTHPtL) 
CALL ECOUN T(l) 




PRINT : J, 
TEM^ TEMP, 
1 C(I) . 



1574 



IL=0 




( PRINT: 

1 

\xc(x+3q)_ 




CALL MSTOfclo) 



IT* I- I 
NOMBERS£l)= 
ICCT+B<rt 




PRIUTjI.KETGV 
KETft*i\lC(ic+«^, 
for i*i>as 




NI»MBERS(P 2 S) 

*IC( 5 l- 75 i 



L 



Box 




CALL 


\0 




HIETO (2$ 
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* Mex"t p3<^€ 
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SUBROUTINE ROUND(L) 



Box 12 
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SUBROUTINE INITLAL(L) 



box 13 
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SUBROUTINE ECOUNT(L) 



Boxes 14 at*) 15 




1= 


1 


J 


L 


Ife40 
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SUBROUTINE HISTO(NGROUPS) 



Box I £> 
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PROGRAM RON 
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APPENDIX V 



DESCRIPTION OF THE FLOW CHARTS 
The explanation given below is to be used with APPENDIX III 
(Program Listings) and APPENDIX IV (Program Flow Charts) . 

PROGRAM MASTER 

BOX 1 

This box writes the title statement, "THIS IS PROGRAM MASTER" 
and as many comment cards (80 spaces per card) as desired. The DO 
loop is written for five comment cards at present, but can be readily 
changed if desired. 

BOX 2 

Reads input parameters for the run. Formats are given in the 
listing in APPENDIX III. 

BOX 3 

Changes lattice dimensions (DC, IY, IZ) into floating point format, 
and computes the physical dimensions of the lattice. These are then 
permanently found in AB (1-3) . It must be kept in mind that IX must 
be ODD . 

BOX 4 

Computes NLINE , NPLANE , MAXLAT , LMAX (the initial LMAX) from 
the formulas given in APPENDIX II. 

BOX 5 

Computes initial lattice positions by use of subroutine SITE. 
Stores the site number and the fact that the site is occupied in the LB 
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array. Stores a -10.0 in each location of the TLAST array. This is 
necessary for collisions occurring in the first time step# otherwise 
these collisions would be ignored. 

BOX 6 

Writes initial lattice positions on Tape 6 in binary for eventual 
use by program RON. Also sets initial or zero values for the indicators 
shown. 

BOXES 7,8, and 9 

Sets initial constant values , sets values of masks , and reads in 
the input parameters for the appropriate Born-Mayer potential, either 
potential 1, 2, or 3. 

BOX 10 

This box includes the regeneration process. Jump Switches 1 & 2 
are raised, leaving Jump Switch 3 down. At the start of the next time 
step, the regeneration information will be written on Tape 7 in binary, 
and the tape will be rewound. To regenerate the run the program must 
be restarted, with all three jump switches in the raised position. After 
the program has started to read Tape 7 the jump switches may be reset. 
When the computer has finished reading Tape 7 , the tape will be re- 
wound. Further regeneration tapes may be made after this point has 
been reached. The main DO loop over the time steps is also started in 
this box. Immediately after the start of a time step the current values 
of TENERGY, N, IENERGY, and ENERGY are written on Tapes 3 & 4. 

The output on Tape 4 is usually transferred to the typewriter to provide 
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a visual check on the validity of the run as it progresses. 

BOX 11 

All decisions concerning continuation or abandonment of the run 
are made in this box. A jump can be made to the output section, the 
normal output dump made, and a return to this box, or abandonment, 
accomplished by the outside control of the jump switches o The run 
will be abandoned, after an output dump, if the most energetic par- 
ticle has an energy less than THERMAL. However, the run will not be 
abandoned under this condition if other particles are to be shot into the 
lattice o 

Additional particles will be shot into the lattice, if the proper jump 
switch is set, at the start of the next time step that is one less than a 
multiple of LULL. For example, if LULL is 3, and the proper switch is 
set in time step 6 (or 7) , an additional particle will enter at the start 
of time step 8. The switch must be left in the set position until after 
the particle has been read into the program. 

The input co-ordinates for the particle are given as - 1.0 for X, 
and the actual point on the front face of the lattice we desire to hit. 

The angles ALPHA & BETA are given, as is the energy. This box then 
calculates the proper co-ordinates and velocity components to achieve 
the desired impact point. Since the new particle is the most energetic 
in the lattice, the Time Step Length is also recalculated at this time. 
BOX 12 

Here, the TSL is recalculated if the number of the time step is an 
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even multiple of JPB , and the first 30 values are placed in the TSLI 
table for output reference . 

BOX 13 

This is the only energy check made during the run, done every time 
step. The total energy (obtained as a sum of the individual energies) 
is compared with the original energy (the total amount shot into the 
lattice) , and if they differ by more than 1%, the MISTAKE (1) error 
indicator is increased by one, 

BOX 14 

The value of ENERGY is recorded in the EMAX table for the first 30 
time steps, the total energy, TENERGY, is recalculated, and the LB 
array is reset and now shows that none of the atoms have been through 
the current time step. 

BOX 15 

The major DO loop on the particles is started here, with the LB 
array used to check certain conditions on each particle. If the atom: 

1) Does not have an energy greater than THERMAL, or 

2) Has left the lattice , or 

3) Has completed this time step, or 

4) Has been through INVAC and has not had a collision since , 
a jump is made to the end of the DO loop, and the next atom is con- 
sidered. 

Therefore, only those atoms are considered that satisfy aU the 
following conditions: 
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1) The atom is a bullet, i.e. has more energy than THERMAL, 

2) The atom is in the lattice , 

3) The atom has not completed this time step, and 

4) The atom has not been through INVAC since its last collision 
If all the preceeding conditions are satisfied, the initial values of the 
counters and lists shown in BOX 15, APPENDIX III are set to zero or other 
appropriate values. M, the fifth Index Register, is then set equal to L, 
the fourth Index Register. L, as such, is not used for the rest of the 
DO loop, until this pass through the DO loop is completed and a new 

L is chosen . 

BOX 16 

All atoms within a mathematical box with sides 3.02 (A) (centered 
at M) are placed on the NUM list here. Atoms are excluded if they have 
left the lattice or have been through the time step. If the list contains 
only one atom (M is placed first on the list) , then a jump is made to BOX 
25 where another value of M is chosen from the LAST list. With other 
atoms in the box, indicators IN3 and NFACE are set to zero. 

BOX 17 

For every atom on the NUM list, the following parameters are com- 
puted with respect to atom M: 

1) The undeviated distance of closest approach between centers , 

called DSTANCE ( ) , 

2) The time of closest approach between centers, also un- 
deviated, called TMIN ( ), 
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3) The time when the two atoms are separated by a distance of 
SPHI, called T ( ). This is the actual time a collision 
would occur. 

The distance is given in Angstroms , with the time in jiffys , relative 
to the start of the current time step. 

Various conditions are now imposed on these quantities to insure the 
validity of these presumed collisions . An atom must pass all the con- 
ditions to insure that a collision will take place. The conditions are: 

1) &V) > 1.0 x 10" 6 , 

2) DSTANCE ( ) SPHI , 

3) T ( ) lies between + TSL of the start of this time step, 

4) T ( ) is less than 0.99999 TSL. This eliminates collisions 
that end too far into the next time step . 

5) The TLAST for both atoms must be less than the absolute 
time of this collision. 

If all of the above conditions are satisfied, IN3 is increased by 
one. If one or more of the conditions is not satisfied, T ( ) is set 

equal to a large positive number (100 TSL serves this purpose in the 
program) . 

It should be noted that DSTANCE (4) , TMIN (4) , and T (4) are the 
parameters associated with the atom that is number 4 on the NUM list, 
i.e. NUM (4), and so on. 

BOX 18 

The minimum time on the T ( ) list is found. All those on the NUM 
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list whose corresponding T's are equal to the minimum time, or within 
(CUTOFF) (TSL) of the minimum time are than placed on the KHIT list. 

Atom M is placed first on the KHIT list. At this point, if INPUT is zero, 
the KHIT list is duplicated on the LAST list. This occurs only when M 
equals L. INPUT is increased and the next box is entered if there is 
more than one atom on the KHIT list, i.e. another atom besides M. If 
there are no other atoms on the KHIT list, a jump is made to BOX 25, 
and M is set equal to the next atom on the LAST list. 

BOX 19 

This box and the next are concerned with the two body interaction 
itself, and the subsequent scaling methods for energy and velocity. 

These scaling methods are used whenever M collides with more than 
one atom. The pre -interaction velocities and energies are stored in 
the SAVE array. The COLD list is set to zero, and J is set to KHIT (JA) . 
JA is the index used for the KHIT list in this box, M then hits the atoms 
in order of their placement on the KHIT list. JA is initially two, and 
increases up to and including NMAX. 

With J selected, LB (J) is checked, and if J is a member of an inter- 
stitial pair, NPAIR is set to 3. The indicator NPAIR is used in BOX 24 
to insure that the appropriate action is taken with respect to J's inter- 
stitial partner after the collision is over. 

Since the atoms on the NUM and KHIT lists do not have the same 
placement, a search is now made to find J on the NUM list. When this 
is done, we have the collision time, T, and the collision parameter, PAR. 
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The two atoms , M & J , are then moved ahead in time to time T . If M 
has already had one or more collisions, it must be moved back in time 
to the start of this time step, and then ahead to the start of this collision. 

The interaction subroutine is then called. Outputs are contained in 
the DEV array as the new positions, velocities, and energies of M & J. 
TLAST is now calculated for J, and J's co-ordinates and velocities 

are set to those of the DEV array. Notice that the TLAST for J does not 

% 

include any portion of the time necessary for the collision. If a portion 
of this time were included, significant collisions would be eliminated. 

The COLD array (1-6) contains the sum of the final position and velocity 
co-ordinates for M as given by the subroutine for all of M's collisions. 
The total energy lost by M (ELOST) then equals the energy gained by 
this J plus that gained by the other J’s on the KHIT list. NHIT shows 
the actual number of atoms that have collided with M . 

BOX 20 

In this box, the scaling, by various means, of the velocities and 
energies of the atoms on the KHIT list is accomplished. 

The temporary quantity T2 is set equal to the previous energy of M 
minus the amount lost by M in all its collisions . 

If T2 is zero or negative*, the energy and velocities of M are set 
to zero. With T2 negative, the energies and velocities of all others 
on the KHIT list must be scaled to conserve energy. All atoms are 
scaled in the same proportion, the absolute direction of their velocity 
vectors is preserved , their energies being scaled by the square of 
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the scale factor used on the velocity components . 



In all cases, the spatial co-ordinates of M after the collision 
must be found. These are set equal to the numerical average of M’s 
positions as computed by the subroutine. This is accomplished by 
simply dividing COLD (1-3) by NHIT . 

For a positive T2, there are two cases , the energy of M before the 
collisions above ETHRESH and below ETHRESH. 

When the original Ej^ ETHRESH and there was only one collision, 
the co-ordinates of M are set to the DEV results, and a jump is made 
to the next box. If there was more than one collision, E m is set equal 
to T2, and the velocities are scaled so that the components will be in 
the same proportion as the sum of the individual resultant velocities 
as computed by the subroutine. This essentially gives M a resultant 
direction equal to the vector sum of the directions M would have as a 
result of its collisions with the rest of the KHIT list. 

When the original E m «^ ETHRESH, we assume that M has been 

oscillating about a site, and its site number is obtained from LB (M) . 

E m is then set equal to T2, and the velocities are scaled so that the 

resultant velocity vector will point toward M’s oscillatory center. The 

scaling factor is (V^/ (AX^ +AY^ + AZ^))^? This is multiplied by 

&X, Ay • an d AZ to obtain the correct components. This scaling auto- 
2 

matically keeps V in proper proportion to E m . 

BOX 21 

TLAST for M is calculated at the start of this box. M's TLAST 
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could not be calculated before this point since M had not yet parti- 
cipated in all its possible collisions. 

This box then insures that all atoms with energies less than 
ETHRESH have site numbers stored in LB, that they are placed on the 
LATER list if their energies are less than THERMAL, or that they are left 
alone if in between the two thresholds . 

If an atom has E^ ETHRESH, its site number must be removed from 
LB (M) , and the bit signitifying an occupied site must be removed from 
LB (site) , unless the atom was once part of an interstitial pair. If it 
was part of an interstitial, then the partner section of LB must be re- 
moved from both LB's. 

If an atom should occupy a site and doesn’t, then a site must be 
found. This is accomplished by going part way through INVAC , INVAC 
is exited by use of the MOVE indicator. If a site is calculated for M, 
and it is occupied, the information pertinent to an interstitial pair is 
recorded in LB. If the site is occupied by an interstitial, M is placed 
on the site anyway, and the fact that a "triple interstitial" was formed 
is recorded by increasing MISTAKE (5). E m is inspected at this point, 
and MISTAKE (3) is increased if E m is negative. The symbolic language 
of this box may appear complicated, but comprehension should be 
facilitated by the flow chart. 

BOX 22 

All other atoms on the KHIT list are now subjected to the same pro- 
cess M underwent in the previous box, except that here a jump is made 
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to Statement 600 to enter INVAC . To preserve the values of I , J, and M, 
and also make M equal J, the values of I, J, and M are stored in MAD, 
NAVY, and LOCATE, respectively. These two boxes might, at some 
future stage of development be combined into one box covering the entire 
KHIT list. However, when this section was developed, the treatment of 
M and the others on the KHIT list was sufficiently different so that this 
combination would have been prohibitively difficult. Most of this dif- 
ficulty has disappeared as the program developed. 

BOX 23 

All atoms on the KHIT list are now backed up, time-wise and space- 
wise, to the beginning of the current time step- Each atom must be 
backed up using its particular T ( ) , since they were not all hit at 

exactly the same time. A search through the NUM list accomplishes 
this, giving the correct value on the T list for the computation. 

BOX 24 

If one of the J’s on the KHIT list was a member of an interstitial 
pair, then NPAIR is three, and a jump is made to Statement 236. The 
relations between atom M (the incoming atom) , atom J (the member of 
the pair struck by M) , and the partner , atom LL (the other member of the 
interstitial pair) , are now considered. 

If J has an energy between the two thresholds (vibrating about its 
site), the information in LB about it must be retained. If the partner, 

LL, has been hit, and has an energy greater than ETHRESH, the MASK1 
portion of both LB’s must be erased. If both have an energy^ ETHRESH, 
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nothing is done to either LB. If, however, J has an energy greater than 
ETHRESH and atom LL has not been hit, we must erase from LB (J) the 
site and partner number, erase from LB (LL) LL's partner number (which 
is J) , and also move atom LL back onto the site. It would be a definite 
error to leave the partner, LL, on its previous split position, so it seems 
more realistic to move it back onto the site. 

BOX 25 

The interactions , scaling , and correction of LB have now been com- 
pleted for this round of collisions. A jump is now made to Statement 200 
in BOX 16. If M has no more collisions this time step, a jump back to 
this point will be made. If M does have further collisions, they are cal- 
culated and carried out. 

Upon entering this box, if IN3 is zero, then no collisions were made 
during the last pass through, and this M has completed the time step. 
The next atom on the LAST list, the INMIN th one, is made atom M. 

INMAX is the limit of this list. INMIN is increased by one, and a jump 
is made to Statement 200. When the LAST list has been exhausted, i.e. 
each one has been through the process until it will have no more colli- 
sions this time step, an exit to the next box is made. 

Any atoms on the LAST list that have been put on the LATER list to 
go through INVAC are not sent back through the collision process since 
they have MASK6 set. 

The name of the atom under primary consideration was changed from 
L to M in BOX 15 to prevent the above process from destroying the value 
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of L, since L is the index on the main DO loop over the atoms . 

BOX 26 

In a DO loop over the LAST list, only those atoms not on the LATER 
list are selected, and then moved to the end of the time step. This fact 
is then recorded in LB. The atom's energy is then compared with 
ENERGY, the largest one becomes the new ENERGY, and INERGY is 
changed if needed. 

BOX 27 

If NZ = 1 , the LATER list is empty, and a jump is made around this 
box. If not, IT1 is set equal to (NZ - 1) , as the upper limit of the DO 
loop to follow. INVAC is then called successively for all of the LATER 
list, writing on Tape 6 in binary the particle number and time step, 
both contained in IU, and the X, Y, and Z co-ordinates. 

BOX 28 

In this DO loop on J, only those atoms on the LAST list that were 
not on the LATER list (and have consequently been through INVAC) are 
considered . 

All three co-ordinates are compared to the lattice dimensions . If 
an atom is outside the physical dimensions of the lattice, a value of 
NFACE is calculated. NFACE and the number of the time step of exit 
are then recorded in LB (M) , the site number and possible interstitial 
partner are erased from LB (M) , and M's number is erased from any 
interstitial partner's LB. The memory of occupancy is erased from the 
site if M is not an interstitial member, and the next atom on the LAST 
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list Is considered. 



This is the end of the major DO loop on L, the next L is chosen, 
or an exit is made to the next box if the DO loop on L is completed, 

BOX 29 

The same procedure is followed here as was followed in BOXES 26 
and 28 above , except that atoms: 

1) That have been through the time step , 

2) Having an energy less than THERMAL, or, 

3) That are outside the lattice 
are not considered. 

This is the end of the major DO loop on N, the time step is completed, 
and either a new time step is started , or an exit to the next box and the 
output section is made. 

BOX 30 

Writes selected lattice and program parameters , and a general lattice 
dump on Tape 8 in binary for eventual use by the SLAVE program. Various 
BCD output is then written on Tape 3 , and if the DO loop on N (the time 
steps) was completed, a (-1) is written in binary on Tape 6 and the pro- 
gram ends. If, however, the DO loop on N was not completed, a jump is 
made back to BOX 11. There, a jump is made back to this box (with the 
resultant (-1) on Tape 6, etc.) if Jump Switch 3 is not set. In the normal 
operational mode, the DO loop is not completed, and the program itself 
terminates computation by a jump from BOX 11 to the start of this box, 
the output is written, the jumps back to BOX 11 and back are made, and 
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the program ends . 

A jump around the BCD section of the output may be made by setting 
Jump Switches 1 & 3. To facilitate resetting of the jump switches, Stop 
Switch 1 may be set. If this is done, the program will stop after writing 
Tape 8 , and again before jumping back to BOX 1 1 . 

Tapes 3 & 4, then, are used for BCD output, and Tapes 6,7, and 8 
are used for binary output over the course of a run (including regeneration) . 
Since Tape 4 output is normally switched to the typewriter, only four out- 
put tapes are needed in the course of a run. The three binary tapes 
(6,7,8) will be rewound at the end of the run, and may be used without 
further action . An END OF FILE mark must be placed on the end of Tape 3 
before printing. 

SUBROUTINE SITE (M) 

BOX 31 

This short subroutine calculates the fixed point co-ordinates of 
the lattice site given as ah input parameter. The co-ordinates are found 
in IB (1-3) after exit from the subroutine. 

SUBROUTINE INVAC (MJ) 

This subroutine, given the number of an atom, either places the atom 
on a vacant lattice site near it, or forms an interstitial pair with one of 
its nearest neighbors . 

BOX 32 

The initial indicator IN 2 is set to zero, and the input atom is 

v 

changed to number MM. The site number occupied by MM, and the fact 



166 



that MM is on the LATER list are now cleared from LB (MM), and the 
fact that MM has gone through INVAC is set in LB (MM) . 

Atom MM's co-ordinates are now inspected to see whether or not MM 
is within one unit of a face of the lattice. If so, the number of this face 
is stored in IFACE (1-3) , and the indicator NIF is increased. The index 
of IFACE is 1 , 2, or 3, depending on whether the face is in the X, Y, or 
Z direction, respectively. 

The co-ordinates of MM are rounded off, and the fixed point co- 
ordinates of the nearest site (or non-site) are found. 

BOX 33 

We must now determine which of the two general cases we have, site 
or non-site, and if IFACE (1-3) is non-zero, which of the 12 special 
cases we are to consider. If the sum of the three fixed point co-ordinates 
is even, then they designate a lattice site, and JM is set to 1 . If the 
sum is odd, then a non-site is designated, and JM is set to zero. The 
sum of the Y and Z co-ordinates is then checked. If odd, JL = 0, and if 
even, JL = -1. If the indicator MOVE is non-zero, a RETURN to the main 
program is executed. 

The formula for IL (the site number) results in an X co-ordinate 
one less or one more than an actual site for the non-site case. The sum 
of JL and JM determine which case. If the sum is -1, then one is added 
to IL. 

At this point, IL is always either: 

1) The site of the rounded off position of MM, or 



167 



2) The site with X co-ordinate one larger than MM's rounded-off 
position. 

The possible special case is now determined. If IFACE (1-3) ^ 0, then 
JB is set equal to the face number with JD = 1, If more than one IFACE 
is non-zero, MISTAKE (2) is increased, since this means a position along 
an edge or comer of the lattice . 

In either case, if IFACE (1-3) 9^0, then JA is non-zero, and JC is 
set equal to 1. If IFACE (1-3) is zero, then JC = JM (either 1 or 0) . 

JB is set equal to the first non-zero IFACE or to zero if IFACE (1-3) 
is zero. JD is then set to 1 if a non-zero IFACE exists, and to zero 
otherwise . 

The case number ( 1 — 1 4) is now established by the formula: 

JW = (JM*6*JD) + 1 + JB + JC . Cases 1 & 2 are the general cases. 

Cases 3-8 are special forms of case 1, with cases 9-14 as special 
forms of case 2. These special cases simply eliminate some of the 
nearest neighbors found in the general case from consideration. 

The lists MOON and NOON are set to zero , and the next box is 
entered . 

BOX 34 

The first section (down to Statement 2451) is a section of computed 
GO TO statements to insure the correct choice of the nearest neighbor 
formulas. The second section places the site numbers of the nearest 
neighbor positions on the MOON list, provided their numbers lie between 
one and MAXLAT. I is the index used for the MOON list, and is always 
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one greater than the number of sites on the list, 

BOX 35 

All of the nearest neighbor numbers found may not have been placed 
on the MOON list (possible negative numbers, or numbers greater than 
MAXLAT) , JPMAX, therefore, is reset to (i-l). If the indicator MOVE is 
non-zero, a RETURN to the main program is executed. 

The LB for each atom whose site is on the MOON list is now in- 
spected. If the LB shows that the site is not occupied (see APPENDIX I) , 
the site number is placed on the NOON list. JP is the index used for the 
NOON list. 

The MOON list, then, is a list of the sites of the nearest neighbors, 
and the NOON list is a list of those nearest neighbor sites that are vacant 
sites, with JK the number of vacancies on the NOON list. 

BOX 36 

If JK = 1, then only one vacancy exists, and atom MM is placed 
there, the site number, IL, is set in LB (MM), and the fact that the site 
is now occupied is stored in LB (IL) . If IN2 is zero, a RETURN is then 
executed, otherwise a jump to BOX 39 (Statement 2800) is executed. 

If JK is zero, there are no vacancies and a prospective interstitial 
partner must be found. If IN2 is zero, a jump to BOX 38 is made to 
accomplish this. If IN2 is non-zero, then MM is the prospective inter- 
stitial partner of the original atom MJ , and a jump is made to BOX 39 to 
form the interstitial. 
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If JK > 1 , then more than one vacancy exists . A jump is made to 
BOX 37, where one of the vacancies is chose, a return to this box is then 
made, and the procedure given for JK = 1 is followed. 

BOX 37 

The (distance) ^ between atom MM and the vacant sites on the NOON 

2 

list is computed and stored in the DSTANCE array. The minimum (distance) 
is then found, and the DSTANCE list is inspected to find the number of 
sites that are at this minimum distance from atom MM. If there is only 
one, IL is set equal to that site, and a jump is made to Statement 2492 in 
BOX 36 . If there is more than one , MM is placed on its own site , if its 
own site is one of those on the list. If not, then MM is placed on the 
first or second site on the NOON list, depending on whether the number 
of the time step is odd or even, respectively. 

The number of the time step is used here as a convenient random 
number generator, since the decision is made on oddness or evenness, 
rather than the actual value of the number itself. 

BOX 38 

This box is entered only from BOX 36, and then on the condition that 
IN2 is zero and that MM (MJ) has found no vacancies. A prospective 
interstitial partner must be chosen from the MOON list. 

The time of atom MM's (MJ's) closest approach to each of the sites 
listed on the MOON list is found and stored in the TMIN list. Any sites 
occupied by interstitial pairs are not considered. 
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The minimum time of closest approach is then found, TMIN (KMIN) , 
which locates the prospective interstitial site. MAX is derived from 
A X, AY, or AZ, whichever is greatest (Z\ X = X co-ordinate of MM 
minus X co-ordinate of the site chosen) , and this becomes the prospective 
axis of the interstitial pair. The atom occupying this site is then found. 
Since it is possible for this atom to be pushed into any vacancy near it, 
and MM to fall into the vacancy thus created, this atom's nearest neigh- 
bors must be inspected. This is done by setting MM equal to the number 
of this atom, setting IN2 to 1, and jumping to Statement 2361 in BOX 32. 
BOX 39 

If the new MM (the original atom is still designated by MJ) does find 
a vacancy, it is placed there in BOX 36, a jump is made to this box, MM's 
now vacant site is named as the prospective site for MJ, a return to 
BOX 36 is made, and MJ is placed on its prospective site. A RETURN to 
the main program is then executed in BOX 36. 

If, however, the new MM did not find a vacant site, then a jump is 
made to this box (Statement 2720), and the interstitial is actually formed. 
The site of the interstitial is entered in each LB, together with the other 
atom’s number. The fact that the site is occupied is then entered in the 
LB of the site, and a RETURN is executed. 

If the new MM equals MJ, then the atom has formed an interstitial 
with itself, and the error counter MISTAKE (4) is increased. The inter- 
stitial is formed so that MJ is on the same side as it was before forma- 
tion of the interstitial pair. 
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SUBROUTINE INTER (M Q, II, S) 

This subroutine, given the atom under consideration (MQ) , the 
target atom (II) , and the impact parameter (S , in meters) calculates the 
positions of the atoms at the end of the collision, their velocities, and 
energies and stores them in the DEV array » 

BOX 40 

Transforms atom MQ's position co-ordinates, velocity components, 
and energy and the angle between the radius vector from atom II to 
atom MQ and MQ's velocity into the target (atom II) frame of reference. 
Sets CPA = SPHI as the initial value for the CPA iteration. 

BOX 41 

Computes CPA by an iterative process. The iteration is terminated 
when two successive values differ by less than 10“^ % , or when any 
value is greater than the preceeding value. 

BOX 42 

Computes the time required for the interaction and the total angle 
of deviation. The integrals are solved by numerical integration using a 
Gauss -Legendre four point quadrature. 

BOX 43 

Computes the composite factors and trigonometric functions that are 
used more than once in the transformation to the laboratory reference 
frame . 

BOX 44 

Computes direction cosines, relative to the laboratory frame of 
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reference, of a line in the plane of the interaction that is perpendicular 
to the line between the two atoms at the beginning of the interaction. 

BOX 45 

Computes the position co-ordinates, velocity components, and 
energies of both atoms and stores them in the DEV array. See APPENDIX I, 
DEV for the actual location. 

PROGRAM SLAVE 

BOX 1 

Reads binary input from Tape 8. "SET SLJ SWITCH 2 UP IF NEEDED" 
is written on Tape 3 (normally transferred to the typewriter) . As many 
comment cards (80 spaces per card) as desired (6 at present) are read 
from Tape 2 and written on the output tape . Reads input parameters for 
the run from Tape 2. Writes on the output tape "SELECTIVE JUMP SWITCH 
NUMBER TWO WAS UP" if Jump Switch 2 is used for the run. 

IHISTO (1-4) and IHISTO (105) are defined, the IHISTO array then 
serves as the format specification for the histographs. The constants 
I2E15, I2E24, F2E24, and I2E30 are set, and the masks for use with the 
LB array are defined . 

BOX 2 

If NTS (called N in PROGRAM MASTER, the number of time steps 
completed) is zero, then NTS is set equal to MNTS; other constants and 
indicators are then set. The times in the TSLI array and the time TSL 
are converted from jiffy s to seconds. Various portions of the LB array 
are then cleared, depending on the settings of the jump switches. 
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BOX 3 



In this box all atoms are placed in one of the following classifica- 
tions: 

1) Atoms occupying their original sites , 

2) Atoms occupying sites other than their original sites , 

3) Atoms that are parts of interstitial pairs , 

4) Atoms whose energy is greater than ETHRESH, and 

5) Atoms outside the lattice. 

Program MASTER has also calculated the face of exit for those atoms that 
have left the lattice . 

The number of atoms and the total energy of those atoms is calculated 
for each classification, as is the number of vacant lattice sites. 

BOX 4 

The total number of atoms and their total energy is calculated. 

This serves as an additional total energy check at the end of the compu- 
tational run. These numbers and various other data (see APPENDIX III) 
are written on the output tape. The number of atoms and their total 
energy is written on the output tape for each classification listed in 
BOX 3. 

BOX 5 

The number of atoms occupying their original sites is written on the 
output tape. An array of these atoms is then written, in order, on the out- 
put tape. Atoms not on their original sites are represented by zeros in 
this array. 
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BOX 6 



The total number of vacant lattice sites and the numbers of these 
sites are written on the output tape. 

BOX 7 

The total number of replacements is written on the output tape. 

The number of the atom, the number of its final site, its penetration, 
and its radial penetration are written on the output tape for each of 
the replacement atoms. Histographs are then compiled and written on the 
output tape for the penetration and the radial penetration distributions 
of the replacement atoms . 

BOX 8 

The total number of interstitial pairs is written on the output tape. 
The number of each atom in a pair, its partner, the site they occupy, the 
penetration of each atom in the pair, and the radial penetration of each 
atom in the pair is written on the output tape. This is written twice for 
each interstitial pair since there are two atoms in the pair. Histographs 
are compiled and written on the output tape for the penetration and the 
radial penetration distributions of all interstitial pair members. 

If Jump Switch 2 was set for the run, the interstitials are determined 
by inspecting the LB array given by program MASTER. If, however, Jump 
Switch 2 was not set, the program disregards the LB array and calculates 
the interstitial pairs in the lattice . 

BOX 9 

The total number of atoms with energy greater than ETHRESH 
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(knock-ons) is written on the output tape . The number of each atom , its 
co-ordinates, velocity components, energy, its LB (in octal), its pene- 
tration , and its radial penetration are written on the output tape for each 
knock-on atom. Histographs are then compiled and written on the out- 
put tape for the penetration , radial penetration , and energy distributions 
of the knock-on atoms . 

BOX 10 

The total number of atoms that exited through a lattice face is 
written on the output tape. The number of the atom, its co-ordinates, 
its velocity components, energy, LB (in octal) , its penetration, radial 
penetration, and time step of exit are written on the output tape for 
each atom exiting through a face . 

The entire box is then repeated for each of the six faces of the 
lattice . 

BOX 11 

If Jump Switch 1 is set, a jump is made to the end of the program, 
otherwise a general lattice dump is written on the output tape. The dump 
consists of the atom number, co-ordinates, velocity components, energy, 
and the LB contents (in octal) for each atom. 

SUBROUTINE ROUND (L) 

BOX 12 

This subroutine computes the site closest to atom L. The atom's 
co-ordinates are rounded off, and the fixed point co-ordinates of the 
atom are stored in the IB array if the co-ordinates represent a site. 
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If the fixed point co-ordinates do not represent a site , the six nearest 
neighbor sites are found. The site that is closest space-wise to the 
atom is chosen as the site for atom L. Once a site is chosen, a RETURN 
to the main program is executed . 

SUBROUTINE INITIAL (L) 

BOX 13 

If atom L is one of the original lattice atoms , this subroutine 
computes its initial fixed point co-ordinates and stores them in the IB 
array. If atom L was shot into the lattice, the point of entry is stored 
in the IB array . A RETURN to the main program is then executed . 

SUBROUTINE ECOUNT (L) 

BOX 14 

The energy of atom L is inspected. This energy will fall between 
two of the values on the KET list, say the sixth and the seventh values 
on the list. Atom L's energy is then greater than the sixth (ith) value 
on the list. The histograph counter IC (I + 50) is then increased by one. 
A RETURN to the main program is then executed. 

SUBROUTINE DEPTH (L) 

BOX 15 

The penetration (DEPTH) is computed as the final X co-ordinate of 
atom L minus the initial X co-ordimte. A value of I is determined so 
that DEPTH lies between A (I -7.5) and A (I -6.5) . Histograph counter 
IC (I) is then increased by one . 
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The radial penetration (PERDIST) is then calculated from L's initial 
and final positions. A value of I is determined so that PERDIST lies be- 
tween A (I -1) and A (I) . Histograph counter IC (I + 30) is then increased 
by one. A RETURN to the main program is then executed. 

SUBROUTINE HISTO (NGROUPS) 

BOX 16 

The largest number and the total of NUMBERS (1 -NGROUpS) are 
determined . The total is written on the output tape . The percentage of 
the total is written on the output tape for each entry on the NUMBERS 
list, followed by a number of X's corresponding to the percentage this 
entry is of the largest number in the table, i.e. the largest number will 
write out 100 X's, and a number half that large will write out 50 X's. 

A RETURN to the main program is then executed. 

PROGRAM RON 

This program writes on the output tape the successive positions of 
every atom in the lattice (or that was shot into the lattice) , if the suc- 
cessive position of an atom is significantly different from the previous 
position. The number of the time step is also written when an atom's 
co-ordinates are written. After writing all the successive positions of 
an atom, the program commences writing the positions of the next atom 
that has moved significantly until all atoms have been considered. 
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